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Abstract 

This  work  is  concerned  with  the  implementation  and  testing  of  several  segregated  multi-fluid  all 
speed  flow  algorithms.  These  algorithms  are  part  of  the  Mass  Conservation  Based  Algorithms 
(MCBA)  group  in  which  the  pressure  correction  equation  is  derived  from  overall  mass  conservation. 
Seven  MCBA  algorithms  (based  on  SIMPLE,  SIMPLEC,  SIMPLEX,  SIMPLEM,  SIMPLEST,  PISO, 
and  PRIME)  are  investigated  within  a  structured  collocated  non-linear  multi-grid  finite-volume 
framework.  The  performance  and  accuracy  of  these  algorithms  are  assessed  by  solving  a  variety  of 
one-  and  two-dimensional  laminar  and  turbulent  two-phase  flow  problems  in  the  subsonic,  transonic, 
and  supersonic  regimes  and  comparing  results  with  published  numerical  and/or  experimental  data. 
The  SG  method  is  used  to  solve  for  the  one-dimensional  test  problems  and  the  effects  of  grid  size  on 
convergence  characteristics  are  analyzed.  However,  solutions  for  the  two-dimensional  problems  are 
generated  for  several  grid  systems  using  the  single  grid  method  (SG),  the  prolongation  grid  method 
(PG),  and  the  full  non-linear  multi-grid  method  (FMG)  and  their  effects  on  convergence  behavior  are 
studied.  The  main  outcomes  of  this  study  are  the  clear  demonstrations  of:  (i)  the  capability  of  all 
MCBA  algorithms  to  deal  with  multi-fluid  flow  situations  when  the  fluidic  continuity  equations  are 
normalized;  (ii)  the  ability  of  the  FMG  method  to  tackle  the  added  non-linearity  of  multi-fluid  flows, 
which  is  manifested  through  the  performance  jump  observed  when  using  the  non-linear  multi-grid  as 
compared  to  the  SG  and  PG  methods;  (iii)  the  extension  of  the  FMG  method  to  predict  turbulent 
multi-fluid  flows;  (iv)  the  capacity  of  the  MCBA  algorithms  to  predict  multi-fluid  flow  at  all  speeds. 
Moreover,  results  displayed  in  terms  of  convergence  history  plots  and  CPU-times,  indicate  that  the 
performances  of  SIMPLE,  SIMPLEC,  and  SIMPLEX  are  very  close.  In  general,  the  performance  of 
SIMPLEST  approaches  that  of  SIMPLE  for  diffusion-dominated  flows  (i.e.  turbulent  flows).  As 
expected,  the  PRIME  algorithm  is  found  to  be  the  most  expensive  due  to  the  explicit  treatment  of  the 
fluidic  momentum  equations.  The  PISO  algorithm  is  generally  more  expensive  than  SIMPLE  and  its 
performance  depends  on  the  type  of  flow  and  solution  method  used.  The  behavior  of  SIMPLEM  is 
consistent  and  in  terms  of  CPU  effort  it  stands  between  PRIME  and  SIMDPLE.  For  all  algorithms,  the 
use  of  the  PG  and  FMG  methods  accelerate  the  convergence  rate.  The  FMG  method  is  by  far  more 
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efficient  than  the  PG  method,  accelerating  the  convergence  rate,  for  the  problems  solved  and  the  grid 
used,  over  the  SG  method  by  a  factor  reaching  a  value  as  high  as  10  and  undoubtedly  verifying  the 
effectiveness  of  multi-grid  methods  in  tackling  the  linearity  of  multi-phase  flows. 
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Nomenclature 

, . .  coefficients  in  the  discretized  equation  for  . 
source  term  in  the  discretized  equation  for  0**^ . 
body  force  per  unit  volume  of  fluid/phase  k. 
coefficient  equals  to  1  /  . 

df  covariant  unit  vector  (i.e.  in  the  direction  of  df ). 

the  D  operator. 

Dp'  the  vector  form  of  the  D  operator, 
drag  function  (Eq.  (4)). 
the  H  operator. 

the  HI  operator  working  on  ((|)®=u*'''y'‘',  or 
the  HP  operator  working  on  ((|)^'‘'=u®y‘',  or  w®). 
HPp[u^*']  the  vector  form  of  HaeHP  operator. 

HIp[u'*']  the  vector  form  of  the  ///operator. 

I**'  inter-phase  momentum  transfer. 

diffusion  flux  of  across  cell  face  T . 

convection  flux  of  across  cell  face  ‘f . 

mass  source  per  unit  volume. 

contravariant  unit  vector  (i.e.  in  the  direction  of  Sf ). 

P  pressure. 

Pr (k)  ^  Pr(k)  jaminar  and  turbulent  Prandtl  number  for  fluid/phase  k. 
heat  generated  per  unit  volume  of  fluid/phase  k. 
general  source  term  of  fluid/phase  k. 
r  volume  fraction  of  fluid/phase  k. 

gas  constant  for  fluid/phase  k. 

Sf  surface  vector, 

t  time. 

temperature  of  fluid/phase  k. 

U^P  interface  flux  velocity  (yP^P^  of  fluid/phase  k. 
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velocity  vector  of  fluid/phase  k. 

^(k)yk)^  velocity  components  of  fluid/phase  k. 

X,  y,z  Cartesian  coordinates. 

||a,  b||  the  maximum  of  a  and  b. 

Greek  Symbols 

density  of  fluid/phase  k. 
diffusion  coefficient  of  fluid/phase  k. 

<E>  dissipation  term  in  energy  equation  of  fluid/phase  k. 

general  scalar  quantity  associated  with  fluid/phase  k. 

Kj-  space  vector  equal  to  (n^  -  yd^  )Sf 

A  p  [0  J  the  A  operator, 

laminar  and  turbulent  viscosity  of  fluid/phase  k. 

Q.  cell  volume. 

thermal  expansion  coefficient  for  phase/fluid  k. 

8t  time  step, 

y  scaling  factor. 

Subscripts 

e,  w, .  refers  to  the  east,  west, . . .  face  of  a  control  volume. 

E,W,..  refers  to  the  East,  West,  . . .  neighbors  of  the  main  grid  point, 

f  refers  to  control  volume  face  f 

P  refers  to  the  P  grid  point. 

Superscripts 

C  refers  to  convection  contribution. 

D  refers  to  diffusion  contribution. 

(k)  refers  to  fluid/phase  k. 

ik)*jik)** ,..  refers  to  first,  second, . . .  updated  value  at  the  current  iteration, 
ik)  o  refers  to  values  of  fluid/phase  k  from  the  previous  iteration. 

(k)'  refers  to  correction  field  of  phase/fluid  k. 

m  refers  to  fluid/phase  m. 

old  refers  to  values  from  the  previous  time  step, 

refers  to  SIMPLEX. 


sx 


Introduction 


The  extensive  developments  that  have  taken  place  in  Computational  Fluid  Dynamics  (CFD)  during 
the  last  three  decades  have  established  this  evolving  technology  as  a  reliable  and  critical  tool  for  the 
simulation  of  a  wide  variety  of  engineering  fluid  flow  processes  (mixing,  solidification,  turbulence, 
.. .).  During  the  past  few  years  several  issues  have  been  addressed.  Concerns  related  to  accuracy  were 
assuaged  through  the  development  of  High  Resolution  (HR)  schemes  [1-13].  Moreover,  better 
solution  algorithms  [14-20],  solvers  [21,22],  and  multi-grid  techniques  [23-26]  have  greatly  reduced 
the  computational  cost  and  made  it  feasible  to  solve  real  life  problems.  Furthermore,  pressure-based 
algorithms  [27-40]  were  extended  to  solve  flow  problems  in  the  various  Reynolds  and  Mach  number 
regimes  using  both  structured  and  unstmctured  grid  methods  [41-46], 

While  high-resolution  schemes,  solvers,  multi-grid  techniques,  structured  and  unstructured  grid 
methods,  etc...  can  be  applied  indiscriminately  to  the  simulation  of  single-fluid  or  multi-fluid  flow, 
nearly  all  developments  in  solution  algorithms  have  been  directed  towards  the  simulation  of  single¬ 
fluid  flow  [24,26,47,48].  In  particular,  many  of  single-fluid  segregated  solution  algorithms  were 
developed  such  as  the  well  known  SIMPLE  [14,15],  the  SIMPLEST  [49],  the  SIMPLEC  [17],  the 
SIMPLEM  [50],  the  PISO  [16],  the  PRIME  [51],  and  the  SIMPLEX  [19]  algorithms,  to  site  a  few. 
Additionally,  several  techniques  have  been  devised  to  improve  the  performance,  facilitate  the 
implementation,  and  extend  the  capability  of  these  algorithms. 

Developments  in  solution  algorithms  for  simulating  multi-fluid  phenomena  have,  on  the  other  hand, 
lagged  behind  that  of  single-fluid  algorithms  due  to  the  much  higher  computational  cost  involved,  to 
numerical  difficulties  that  had  to  be  first  addressed  in  the  simulation  of  single-fluid  flow,  and  to  the 
increase  in  algorithmic  complexity.  While  the  major  difficulty  in  the  simulation  of  single-fluid  flow 
stems  from  the  coupling  between  the  momentum  and  continuity  equations,  in  the  simulation  of  multi¬ 
fluid  flow  phenomena,  this  problem  is  further  complicated  by  the  fact  that  there  are  as  many  sets  of 
continuity  and  momentum  equations  as  there  are  fluids,  that  they  are  all  coupled  together  in  various 
ways  (interchange  momentum  by  inter-phase  mass  and  momentum  transfer,  etc.)  and  that  the  fluids 
share  space  (the  volume  fractions  sum  to  unity,  but  are  not  known  in  advance). 
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Despite  these  complexities,  successful  segregated  pressure-based  solution  algorithms  have  been 
devised.  The  IPSA  variants  devised  by  the  Spalding  Group  at  Imperial  College  [52-54]  and  the  set  of 
algorithms  devised  by  the  Los  Alamos  Scientific  Laboratory  (LASL)  group  [55-57]  are  examples  of 
multiphase  algorithms.  However,  in  contrast  with  the  widespread  information  available  on  single¬ 
fluid  solution  algorithms,  little  information  is  available  about  multi-fluid  solution  algorithms,  a  fact 
that  has  confined  their  use  to  a  small  community,  slowed  their  development,  and  isolated  them  from 
newer  developments  in  single-fluid  flow  algorithms  (PWIM,  all  speed  flows,...). 

Recently,  Darwish  et  al.  [58,59]  extended  the  segregated  single-fluid  flow  algorithms  [20]  to  predict 
multi-fluid  flow  at  all  speeds  and  derived  these  algorithms  using  a  unified,  compact,  and  easy  to 
understand  notation.  Furthermore,  it  was  shown  in  [58,59]  that  the  pressure  correction  equation  can  be 
derived  either  by  using  the  geometric  conservation  equation  or  the  overall  mass  conservation 
equation.  Depending  on  which  equation  is  used,  the  segregated  pressure-based  multi-fluid  flow 
algorithms  were  classified  respectively  as  either  the  Geometric  Conservation  Based  family  of 
Algorithms  (MCBA)  or  the  Mass  Conservation  Based  family  of  Algorithms  (GCBA). 

The  objective  of  the  present  work  is  to  test  and  compare  the  different  algorithms  in  the  MCBA  group 
and  to  assess  their  performance  by  solving  various  one-  and  two-dimensional  laminar  and  turbulent 
two-phase  problems  in  the  subsonic,  transonic,  and  supersonic  regimes.  In  the  implementation  of 
these  algorithms  care  is  taken  to  use  a  high  accuracy  discretization  within  the  framework  of  a  non¬ 
linear  full  multi-grid  solution  procedure  [60-66]. 

In  what  follows,  the  governing  equations  are  first  introduced,  followed  by  a  brief  description  of  the 
discretization  procedure.  Then  the  MCBA  algorithms  are  presented,  their  capabilities  to  predict  multi¬ 
fluid  flow  phenomena  at  all  speeds  demonstrated,  and  their  performance  characteristics  (in  terms  of 
convergence  history  and  speed)  assessed.  For  that  purpose,  a  total  of  twelve  laminar  and  turbulent 
incompressible  and  compressible  problems  encompassing  dilute  and  dense  gas-solid,  and  bubbly 
flows  in  the  subsonic,  transonic  and  supersonic  regimes  are  solved.  In  addition,  the  performance  of 
these  algorithms  on  a  single  grid  and  within  a  multi-grid  environment  is  compared. 
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The  Governing  Equations 

In  multi-fluid  flow  the  various  fluids/phases  coexist  with  different  concentrations  at  different 
locations  in  the  flow  domain  and  move  with  unequal  velocities.  Thus,  the  equations  governing  multi¬ 
fluid  flows  are  the  conservation  laws  of  mass,  momentum,  and  energy  for  each  individual  fluid  in 
addition  to  a  set  of  auxiliary  relations. 


Conservation  of  mass 


The  volume  fraction  which  is  the  portion  of  volumetric  space  occupied  by  the  k*  flmd  (Q^'^VQ) 
along  with  the  k'"'  fluid  density,  p®,  and  velocity,  u®,  in  order  to  satisfy  the  mass-conservation 
principle,  have  to  obey  the  differential  equation: 


at  ^  ^ 

Mass  sources  are  often  non-zero,  as  when  one  fluid  is  transformed  to  another  fluid.  However, 
summation  over  all  fluids  leads  to  the  following  “global  mass-conservation”  eqxxation: 


at 


=  0 


(2) 


V  ■  / 

The  zero  on  the  right-hand  side  signifies  that  the  sum  of  mass  sources  (generation  and  loss)  is  zero. 


Conservation  of  momentum 

Denoting  the  velocity  of  the  k'*’  phase  by  then  the  fluidic  momentum  equation  (for  the  k'*'  phase) 
can  be  written  as: 

(k)  {k)\ 

Here  P  stands  for  the  pressure,  which  is  assumed  to  be  shared  by  all  fluids,  is  the  body  force  per 
unit  volume  of  phase  (k),  and  |xf  We  the  laminar  and  turbulent  viscosity,  and  represents 
the  interfacial  forces  per  unit  volume  due  to  drag,  virtual  mass  effects,  lift,  etc.  Because  of  the  many 
suggested  closures  for  the  various  interfacial  sources,  the  ones  adopted  in  the  present  study  will  be 
presented  as  needed. 
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Conservation  of  energy 

Let  be  the  temperature  of  the  k'*'  phase,  the  fluidic  energy  equation  for  the  k***  fluid  is  then  given 


by: 


V.(r®p<'' ¥■' >!'''>)=  V. 


Pr 


,(k) 


Pr; 


.(k) 


VT' 


,(k) 


J 


,.(k) 


,  ^ + v.(pu<'‘>)-  pv.(u<''>)  +  + ^''4 +44 

L  at  J  J  Cp 


(4) 


where  is  the  viscous  dissipation  function,  the  thermal  expansion  coefficient  which  is  equal 
to  for  an  ideal  gas,  Pr^^^  and  Prf^^are  the  laminar  and  turbulent  Prandtl  numbers,  and  ^is 
the  interfacial  energy  transfer  to  phase  (k). 

Turbulence  Modeling 

The  effect  of  turbulence  on  interfacial  mass,  momentum,  and  energy  transfer  is  difficult  to  model  and 
is  still  an  active  area  of  research.  Similar  to  single-fluid  flow,  researchers  have  advertised  several 
flow-dependent  models  to  describe  turbulence.  These  models  vary  in  complexity  from  simple 
algebraic  [67]  models  to  state-of-the-art  Reynolds-stress  [68]  models.  However,  the  widely  used 
multi-phase  turbulence  model,  presented  next,  is  the  two-equation  k-£  model  [69].  Without  going  into 
details,  the  fluidic  conservation  equations  governing  the  turbulence  kinetic  energy  (k)  and  turbulence 
dissipation  rate  (e)  for  the  k^^  fluid  are  given  by: 


^  y^(^(k)p(k)^(k)^(k))^ 


at 

3^,(k)p(k)g(k)) 

at 


^  ,|(k)  ^ 

r«^Vk» 

V  ^k  J 


(k) 

k 


+ 


V.(r(k)p(k)u(k)e(k))^  V. 


,.(k) 


^.(k) 


r(k) 


(5) 


(6) 


r,(k) 


-j^  (ci eG®  -  C2,e<‘=>)  +  If  ^ 

where  and  represent  the  interfacial  turbulence  terms.  The  turbulent  viscosity  is  calculated  as: 


(k)_c 


+  e(k) 


(V) 
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For  two-phase  flows,  several  extensions  of  the  k-s  model  that  are  based  on  calculating  the  turbulent 
viscosity  by  solving  the  k  and  B  equations  for  the  carrier  or  continuous  phase  only  have  been  proposed 
in  the  literature  [70-76].  In  a  recent  article,  Cokljat  and  Ivanov  [69]  presented  a  phase  coupled  k-£ 
turbulence  model,  intended  for  the  cases  where  a  non-dilute  secondary  phase  is  present,  in  which  the 
k-B  transport  equations  for  all  phases  are  solved.  Since  the  method  is  still  not  well  developed,  the  first 
approach  in  which  only  the  k  and  B  equations  for  the  carrier  phase  are  solved  is  adopted  in  this  work. 
Details  regarding  the  specific  model  used  will  be  presented  as  needed. 

The  General  Fluidic  Scalar  Equation 

A  review  of  the  above  differential  equations  reveals  that  they  are  similar  in  structure.  If  a  typical 
representative  variable  associated  with  phase  (k)  is  denoted  by  the  general  fluidic  differential 
equation  may  be  written  as: 

where  the  expression  for  and  can  be  deduced  from  the  parent  equations. 

Auxiliary  Relations 

The  above  set  of  differential  equations  has  to  be  solved  in  conjunction  with  constraints  on  certain 
variables  represented  by  algebraic  relations.  These  auxiliary  relations  include  the  equations  of  state, 
the  geometric  conservation  equation,  and  the  interfacial  mass,  momentum,  energy,  and  turbulence 
energy  transfers. 

Physically,  the  geometric  conservation  equation  is  a  statement  indicating  that  the  sum  of  volumes 
occupied  by  the  different  fluids  within  a  cell  is  equal  to  the  volume  of  the  cell  containing  the  fluids. 

(9) 

k 

For  a  compressible  multi-fluid  flow,  auxiliary  equations  of  state  relating  density  to  pressure  and 
temperature  are  needed.  For  the  phase,  such  an  equation  can  be  written  as: 

p(k)^p(k)(pTa)) 


(10) 
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Several  models  have  been  developed  for  computing  the  interfacial  mass,  momentum,  energy,  and 
turbulence  energy  transfers  terms.  The  closures  used  in  this  work  will  be  detailed  whenever  they  arise 
while  solving  problems. 

In  order  to  present  a  complete  mathematical  problem,  thermodynamic  relations  may  be  needed  and 
initial  and  boundary  conditions  should  supplement  the  above  equations. 

Discretization  Procedure 


In  the  previous  sections  the  differential  equations  governing  multi-fluid  flow  phenomena  were 
outlined  as  well  as  the  associated  auxiliary  relations.  The  task  now  is  to  present  the  Finite  Voliune- 
based  numerical  solution  algorithm  for  solving  these  equations. 


Discretization  of  the  general  fluidic  conservation  equation 

The  general  conservation  equation  (8)  is  integrated  over  a  finite  volume  to  yield: 


ff  fV.(r<">p*u»>0«)iQ 

n  *  a  (11) 


=  ff  V.(r<‘'^r<'‘V(t)<''))dQ  +  Jj  r<‘'^Q^''Ma 

n  o 

Where  Q.  is  the  volume  of  the  control  cell.  Using  the  divergence  theorem  to  transform  the  volume 
integral  into  a  surface  integral  and  then  replacing  the  surface  integral  by  a  summation  of  the  fluxes 
over  the  sides  of  the  control  volume,  equation  (1 1)  is  transformed  to: 

X(jr+J‘i*)=r”Q"'ii  (12) 

nb=e,w,  n,s,t,b 

where  and  are  the  diffusive  and  convective  fluxes,  respectively.  The  discretized  form  of 
the  diffusive  flux  along  an  east  face  is  given  by: 


where  the  over  bar  denotes  a  value  obtained  by  interpolation,  is  a  space  vector  defined  as, 


(13) 


Ke  =  +  <i  + 

and  Y  is  a  scaling  factor  given  by  [77], 


(14) 
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“ede  S,.d, 

and  n^and  df  are  the  contravariant  and  covariant  unit  vectors  respectively.  The  discretized 
convective  flux  along  an  east  face  is  given  by: 


,(k)C 


,(k),(k) 


(16) 


where  Seand  Ce  are  the  surface  vector  and  convection  flux  coefficient  at  cell  face  'e',  respectively.  As 
can  be  seen  from  Eq.  (16),  the  accuracy  of  the  control  volume  solution  for  the  convective  scalar  flux 
depends  on  the  proper  estimation  of  the  face  value  (|)e  as  a  function  of  the  neighboring  (|)  nodes  values. 
Using  some  assumed  interpolation  profile,  (|)e  can  be  explicitly  formulated  in  terms  of  the  nodal  values 
by  a  functional  relationship  of  the  form: 


=  O’) 

where  denotes  the  (j)  values  at  the  neighboring  nodes.  The  interpolation  profile  may  be  one¬ 
dimensional  or  multi-dimensional  of  low  or  high  order  of  accuracy.  The  higher  the  order  of  the  profile 
is,  the  lower  numerical  diffusion  will  be.  However,  the  order  of  the  profile  and  its  dimensionality  do 
not  eliminate  numerical  dispersion.  To  minimize  this  error,  limiters  on  the  convective  flux  should  be 
imposed  to  enforce  monotonicity.  The  flux  limiter  denoted  by  the  Convective  Boundedness  Criterion 
(CBC)  [2]  is  adopted  here  and  applied  within  the  context  of  the  Normalized  Variable  and  Space 
Formulation  methodology  (NVSF)  [8]  on  a  quadratic  upwind  biased  profile,  which  is  third  order 
accurate,  to  produce  the  SMART  [2]  convective  scheme.  Unless  otherwise  stated,  the  High  Resolution 
(HR)  SMART  scheme  is  used  in  solving  all  problems  presented  here.  For  more  details  regarding  the 
NVSF  methodology  the  reader  is  referred  to  [8]. 

After  substituting  the  face  values  by  their  functional  relationship  relating  to  the  node  values  of  cj),  Eq. 
(12)  is  transformed  after  some  algebraic  manipulations  into  the  following  discretized  equation: 

Ai«tf>=lA®C+B'«  (18) 

NB 

where  the  coefficients  and  A^  depend  on  the  selected  scheme  and  is  the  source  term  of 
the  discretized  equation  .  In  compact  form,  the  above  equation  can  be  written  as 
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(!><'''  =  Hp[(t)<''^> 
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Discretization  of  the  fluidic  momentum  equation 

The  discretization  procedure  for  the  momentum  equation  yields  an  algebraic  equation  of  the  form: 

A'J’ur  =  I  A»>««  +  B«  -  r;‘>n,V,(P)+  Q,  Xg<‘">«’  -  u«)  (20) 

NB(P)  m^k 

In  the  above  equation,  the  inter-phase  term  is  written  out  explicitly  to  show  the  strong  coupling 
among  the  momentum  equations  of  the  different  fluids.  This  is  different  from  the  spatial  coupling  that 
exists  among  the  neighboring  velocities  of  the  same  fluid.  To  improve  the  overall  convergence  and 
robustness  of  the  algorithm,  the  inter-phase  source  term  can  be  linearized  and  the  linear  part  treated 
implicitly  so  that: 

A',‘’u'‘>  =  X  aSu®  +  B«  -  r®n,  V,(P)  +  £1,  Xg'‘«u“  (21) 


Where  now 

A';>^A®  +  Q,X(gl,'"')  (22) 

m9£k 

For  later  reference,  the  discretized  form  of  the  momentum  equation  is  written  in  the  following  two 
forms: 

u‘''^  =  HPp[V'‘^]-r®D|!'Vp(P)  (23) 

where  the  body  force  and  inter-phase  terms  are  absorbed  in  the  Bp*'  'source  term  within  the  HP p  [u^  '  ] 

term,  or  as 

u®=Hi,[«“]+Dri(g“"’iin  <24) 

where  the  body  force  and  pressure  gradient  terms  are  absorbed  in  the  Bp'''  source  term  within  the 
Hip  [u^''']  term. 

Discretization  of  the  fluidic  mass  conservation  equation 

The  fluidic  mass-conservation  equation  (Eq.  (1))  can  be  viewed  as  a  fluidic  volume  fraction  equation: 

r‘'''=Hp[r<''']  (25) 

or  as  a  fluidic  continuity  equation  to  be  used  in  deriving  the  pressure  correction  equation: 
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Q+Ap[r^^>p^'‘>u«.S]=r<^^#^ 


where  the  ?  operator  represents  the  following  operation: 


Ap[0]=  £0f 

f=NB{P) 

Discretization  of  the  fluidic  energy  equation 


(26) 


(27) 


The  discretization  of  the  energy  equation  follows  that  of  the  general  fluidic  scalar  equation.  The  only 
difference  is  the  one  pertaining  to  the  discretization  of  the  additional  source  terms.  Since  a  control 
volume  approach  is  followed,  the  integral  of  these  sources  over  the  control  volume  appears  in  the 
discretized  equation.  By  using  the  divergence  theorem,  the  volume  integral  is  transformed  into  a 
surface  integral  and  the  resultant  discretized  expressions  evaluated  explicitly. 


The  Non-linear  Multi-Grid  Strategy 

Standard  relaxation  techniques  reduce  those  Fourier  components  of  the  error  whose  wavelengths  are 
comparable  with  the  mesh  size,  and  then  performance  degrades  when  dealing  with  the  low-frequency 
components  of  the  error.  The  underlying  idea  of  the  multi-grid  strategy  is  to  use  progressively  coarser 
grids  on  which  the  low-frequency  error  components  on  the  finest  grid  appear  as  high-frequency 
Fourier  mode  for  which  the  relaxation  scheme  works  efficiently.  In  the  present  work,  this  strategy  is 
adopted  to  accelerate  convergence  and  thereby  reduce  the  overall  computational  cost. 

The  full  multi-grid  cycle  adopted  in  this  work  can  be  summarized  as  follows.  The  coarser  grids  are 
constructed  by  agglomerating  four  fine  grid  cells  in  a  square  topology.  The  algorithm  starts  at  the 
coarsest  level  and  the  solved  fields  are  interpolated  onto  the  finer  mesh  and  used  as  initial  guess.  Few 
iterations  are  performed  on  the  fine  mesh  before  the  solution  is  restricted  back  to  the  coarser  mesh  to 
be  used  as  an  initial  guess  for  the  solution  at  this  level.  A  forcing  term  is  added  to  the  discrete 
conservation  equations  on  the  coarser  grid  to  ensure  that  the  equation  at  the  coarse  grid  represents  the 
fine  grid  equation.  After  performing  a  number  of  iterations  on  the  coarse  mesh,  the  solution  is 
prolongated  onto  the  finer  mesh  and  a  number  of  iterations  are  completed.  This  process  continues 
until  convergence.  The  solution  is  then  interpolated  into  a  finer  mesh  and  the  process  repeated  imtil 
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convergence  is  reached  on  the  desired  finest  mesh.  This  strategy  has  been  applied  to  both 
incompressible  and  compressible  multi-fluid  flows. 

The  Mass  Conservation  Based  Algorithms 

The  number  of  equations  describing  an  n-fluid  flow  situation  are:  n  fluidic  momentum  equations,  n 
fluidic  volume  fraction  (or  mass  conservation)  equations,  a  geometric  conservation  equation,  and  for 
the  case  of  a  compressible  flow  an  additional  n  auxiliary  pressure-density  relations.  Moreover,  the 
variables  involved  are  the  n  fluidic  velocity  vectors,  the  n  fluidic  volume  fractions,  the  pressure  field, 
and  for  a  compressible  flow  an  additional  n  unknown  fluidic  density  fields.  In  all  MCBA  algorithms, 
the  n  momentum  equations  are  used  to  calculate  the  n  velocity  fields,  n-1  volume  fraction  (mass 
conservation)  equations  are  used  to  calculate  n-1  volume  fraction  fields,  and  the  last  volume  fraction 
field  calculated  using  the  geometric  conservation  equation 

(28) 

The  remaining  volume  fraction  equation  can  be  used  to  calculate  the  pressure  field.  However,  instead 
of  using  this  last  volume  fraction  equation,  the  global  conservation  equation  can  be  employed,  i.e.  the 
sum  of  the  individual  mass  conservation  equations,  to  derive  a  pressure  correction  equation.  The 
sequence  of  events  in  the  MCBA  is  as  follows: 

1.  Solve  the  fluidic  momentum  equations  for  velocities. 

2.  Solve  the  pressure  correction  equation  based  on  global  mass  conservation. 

3.  Correct  velocities,  densities,  and  pressure. 

4.  Solve  the  fluidic  mass  conservation  equations  for  volume  fractions. 

5.  Solve  fluidic  the  energy  equations. 

6.  Return  to  the  first  step  and  repeat  until  convergence. 

Solving  for  Velocities 

In  this  first  step,  the  fluidic  momentum  equations  are  solved  to  find  based  on  guessed  or 

previously  calculated  volume  fraction  and  pressure  fields: 


yi  j  /  c-x^uocu  jt  ■  r  jwci/tty«  vviift  i  uii  iviutti-yjt  tu  uttun  jLft  ivium-i  tutu  i  iuvv  ui /tlh  u^c-c^u^ 


=  fflp  [u®  ]+  (29) 

From  Eq.  (29)  it  is  clear  that  the  HIplu^''^]  term  couples  Up'‘^  to  the  neighbouring  phase  (k) 

velocities  (geometric  or  spatial  coupling  while  hie  Dp''^^g^''"’^Up“^  term  couples  Up*'^  to  the 

velocity  of  all  other  phases  at  grid  point  P  {inter-phase  coupling).  Therefore,  the  rate  of  convergence 
of  the  iterative  solution  procedure  used  to  solve  the  above  system  will  greatly  depend  on  its  capability 
to  resolve  both  types  of  coupling.  The  spatial  coupling  presents  no  problem  to  the  well-established 
iterative  techniques  since  it  couples  velocities  of  the  same  phase.  The  inter-phase  coupling  is  however 
problematic  since  it  relates  velocities  of  different  phases.  An  explicit  evaluation  of  this  term  slows  the 
convergence  rate  considerably  especially  when  the  inter-fluid  momentum  transfer  terms,  represented 
by  g^””^  ,  are  large.  To  accelerate  convergence,  the  Partial  Elimination  Algorithm  (PEA)[78]  or  the 
Simultaneous  solution  of  Non-linearly  Coupled  Equations  technique  (SINCE)  [79]  is  used. 


Solving  for  Pressure  Correction 

To  derive  the  pressure-correction  equation,  the  mass  conservation  equations  of  the  various  phases  are 
added  to  yield  the  global  mass  conservation  equation  given  by: 


Vp  pp  ;  pp ;  ^  ^  u^''>.s) 


(30) 


In  the  predictor  stage  a  guessed  or  an  estimated  pressure  field  from  the  previous  iteration,  denoted  by 

P  is  substituted  into  the  momentum  equations.  The  resulting  velocity  fields  denoted  by  u  which 
now  satisfy  the  momentum  equations  will  not,  in  general,  satisfy  the  mass  conservation  equations. 
Thus,  corrections  are  needed  in  order  to  yield  velocity  and  pressure  fields  that  satisfy  both  equations. 
Denoting  the  corrections  for  pressure,  velocity,  and  density  by  P',  u^’‘^ ,  and  respectively,  the 
corrected  fields  are  written  as: 


P  =F  +P',u^'‘>  =  p®°  +  p<''''  (31) 

Hence  the  equations  solved  in  the  predictor  stage  are: 

u^’‘>*  =  HPp[«®*]-r*“D|.XP° 

While  the  final  solutions  satisfy 


(32) 


n  i  r  K>.3JUf  c-uu>icu  i  ‘  trine,  r  uiurrie  ivietnuu  vviui  i  ttn  ivium-'^jr  lu  JT.eeeiet  uttuti  jur  iviuiti-i  tutu  i  tuvv  ui  Jiti  ijyeeu^ 


±  / 


=  HPp[u<''’]  VpP  (33) 

Subtracting  the  two  equation  sets  ((33)  and  (32))  from  each  other  yields  the  following  equation 

involving  the  correction  terms: 

=  HPp[u*'‘’']  -  r^"^"D|,‘'VpP'  (34) 

Moreover,  the  new  density  and  velocity  fields,  and  will  satisfy  the  overall  mass  conservation 

equation  if: 


6t 


Old 


■a+Ap[r<^V'V'^\s] 


=  0 


Linearizing  the  (p^^u®)  term,  one  gets 


(35) 


+  p<'‘' +  p'*'^  (36) 

Substituting  equations  (36)  and  (34)  into  equation  (35),  rearranging,  and  replacing  density  correction 

by  pressure  correction,  the  final  form  of  the  pressure-correction  equation  is  written  as: 


Ap  [r<'‘^°U“‘>*C5p''^P']-  Ap  D^‘' Vp')s  J 


(k)0  (kr_  /  (k)  (k) 

Tp  Pp  \Lp 

5t 

[+ Ap (hp[u®'])s]+  Ap p^’^’u^'.s] 

The  second  order  correction  term  p^''^  u^'‘^  is  usually  neglected.  Moreover,  if  the  HP  [ii^^  ]  term  in 

the  above  equation  is  retained,  there  will  result  a  pressure  correction  equation  relating  the  pressure 
correction  value  at  a  point  to  all  values  in  the  domain.  To  facilitate  implementation  and  reduce  cost, 
simplifying  assumptions  related  to  this  term  have  been  introduced.  Depending  on  these  assumptions, 
different  algorithms  are  obtained.  A  summary  of  the  various  GCBA  algorithms  used  in  this  work  is 
given  next. 


nOM 
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(k)* 


(37) 


The  MCBA’s  in  Symbolic  Form 

In  what  follow,  the  multi-fluid  versions  of  the  single  fluid  SIMPLE  [14,15],  SIMPLEC  [17], 
SIMPLEST  [49],  SIMPLEX  [19],  SIMPLEM  [50],  PISO  [16],  and  PRIME  [51]  segregated  algorithms 
are  given  using  a  symbolic  format.  In  the  notation  used,  the  superscripts  “old”  and  “o”  denote  values 
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from  the  previous  time  step  and  values  from  the  previous  iteration,  respectively.  Moreover,  the 
superscripts  *,  **,  ***,  and  ****  represent  the  first,  second,  third,  and  fourth  updated  values  at  the 
current  iteration,  respectively. 

The  MCBA  following  SIMPLE  (MCBA-SIMPLE):  Symbolic  Form 

Predictor: 

u«*  =  (38) 

Corrector: 

(u^*^**  =  =  P  =  P<^^”  d-p^'^O  (39) 

=  HP^[u^''^*  ]  -  +P')  (40) 

"|p®'=c^*y 

Condition: 


=  -Si  Si  - -  ^  i 

*  [+A^p*>'p‘**°  (hp[u‘''’'])s]+ p‘*''u^*’'.s]J 


Approximation: 

Neglect: 


HP  P^'],  App^’ p'^^'u^^^S]  =» 


Approximate  Equation: 
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A  Global  MCBA-SIMPLE  Iteration 


Solve  implicitly  for  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P  and 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. 
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The  MCBA  following  SIMPLEC  (MCBA-SIMPLEC):  Symbolic  Form 

Predictor: 

Corrector: 

+u^">‘  =P°  +P^'0 

:.nf"  =HP^[u^''^]-r^^^°D^*V;,P’  =  HP^[u®*  +P') 

.-.  uf'  =  HPp[u‘"^']- 

Subtracting  HPp[l]u®  from  both  sides,  one  gets 
=  HPp[u‘^^']-HPp[l]u^'^' 

'ay  P' 
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[pm-,C';>P 

Condition: 
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Approximate  Equation: 
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A  Global  MCBA-SIMPLEC  Iteration 
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'  Solve  implicitly  for  u(k),  using  the  oid  pressure  and  desity  fieids. 

^Caicuatethe  fields. 

•  Solve  the  pressure  correction  equation. 

>  Correct  P  and 

>  Soive  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fieids. 

'  Return  to  the  first  step  and  iterate  until  convergence. _ 


The  MCBA  following  PRIME  (MCBA-PRIME):  Symbolic  Form 


Predictor: 


=  HP4u®“] 

Corrector: 

/.  nT  =  HP.Cu^'^’**]  =  HP,[u<‘=>*  +  u'*>']  -  +P') 

Uif'  =  HPp[u^"'*  -  +  HP^[u<*^']  -  r^’^^°Y}fVpF 

=  Cj*>P' 

Condition: 
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Approximation: 


Neglect:  HP[u'*>-  -u""],  HP[««]  A,P*>>'‘'ll'‘’.s] 


Approximate  Equation: 
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A  Global  MCBA-PRIME  Iteration 


'  Solve  explicitly  for  vr  \  using  the  oid  pressure  and  density  fields. 

■  Calculate  the  fields. 

<  Solve  the  pressure  correction  equation. 

•  Correct  u^’‘\  P  and  p^l 

'  Solve  implicitly  for 

'  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

<  Return  to  the  first  step  and  iterate  until  convergence. _ 
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The  MCBA  following  SIMPLEST  (MCBA-SIMPLEST):  Symbolic  Form 


Predictor: 


Corrector: 

(u«"  +  u'^>>*  =P“  +P',p^^’‘  =  p‘^^”  +p<^0 

=  HPp'’ [u<''^*]+  HPp°[u^'')  ]+  HPp^[u^''^ ]+  HPp''  ] 

-  t<'')“D(p'‘)VpP“  -  rW-Dj^VpP' 

\uf'  =  HP^[u^'n+  HPp^[u<^^’  -  rW°D<;VpP' 


pW'^Ci*)p- 

Condition: 


Ir  LOf 


=-S: 


rypy-irypy)  , 


& 


[+A,  (hP^[u'*>*  -  ]+  HP[u(*>'])s]+  A, p*)  p^*>'u‘*^'.S]J 

Approximation: 

Neglect: 

HPp^'^'],  =>  uy  ^-y^^jyyvpF 

Approximate  Equation: 

r«“c<;>p; + A^p'^^°c<;^t/<*’*p']-  a^  D('^)vp')st 

^  I5r  4 

[  /k>^(ky  _  (M)  (k)  'f‘‘‘  1 

=  _y  J  IL_££ — Pp  )  Q + A  p  1  }■ 

rj  &  ^  1 

A  Global  MCBA-SIMPLEST  Iteration 


•  Solve  for  treating  diffusion  implicitly  and  convection  explicitly. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P  and  p^l 

•  Solve  implicitly  for  y\ 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 


(60) 

(61) 

(62) 

(63) 


(64) 


(65) 


(66) 


The  MCBA  following  PISO  (MCBA-PISO):  Symbolic  Form 

First  Predictor: 


u®*  =  HPp[u^^  ]  -  y^^DfVpF 


(67) 
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First  Corrector: 

=P°  +P',p'^'’  =  p*^'°  +p'^0 

.-.  =  HPp[u<^»**]  -  -  r'*>“D'/V^(P”  +  F) 

\uf’  =  HP;,[u<*>']-  r<*>°DWVpP' 

•••jpW'^Ci^y 

Condition: 

X  |-  A,  P*'”  D'*  VP')s]j 

t  Lor  -' 


=-i^ 


4«-pi,‘--(r»>pi,‘>r. 


» 


-Q+A,p‘>y‘’‘f^“’'] 


+A,  p’’  p‘*>’  ^[u'‘>])sp  A,[r‘*’'p“''n®'.s]J 

Approximation: 

Neglect: 

HP  [u®'  ],  A,P‘>-p'*V‘’.s]  ^  «f  =  -r'*'-D»'V,P' 

Approximate  Equation: 

2;  j^rW»C</)p;  +  A^  p*) "pW«p*)  D'^>VP')st 

K  [^St  "J 


*  [ 


f//’-p®-(//»p®r. 


•  a+A,pV*'-c/<‘'']l 

J 

Second  Corrector: 

(u«',P",p'‘'' )  (u'*'*"  =  u®"’  +  u“’',  p"  =  p'  +P",p'‘>"  =p“’'  +  p“>') 

.■.  u<;>~"  =  HP;‘tu<‘'”"  ]  -  r">-D«>"V,(p'  +  P") 


.(*)“* 


uf  =Hpr[u^ 

p(k)'  _  Q(lc)p" 


.uW**  +  „«']_rW”D(''^*VpP" 


Condition: 


a 


.(*) 


&  ^ 


‘C<;'P"+A,[r'*>'C<;W‘>'“P"]-A,p>'p'*’'(''‘*’''’'‘’''’''')®} 


=-n 


Old 


5t 


+A^  P*V<*^‘(hP**  +  u<*''])s]+  A,  p*' p<*'"u'^^s]J 


(68) 

(69) 

(70) 


(71) 


(72) 


(73) 


(74) 

(75) 

(76) 

(77) 


(78) 
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Approximation: 

Neglect: 


(kY 


HP"  P®"u^'^'.S] 

Approximate  Equation: 

^  ^  ppJ_q+ 


=-S| 

*  [ 


St 


J 


A  Global  MCBA-PISO  Iteration 


(79) 


(80) 


•  Solve  implicitly  for  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,and 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equation  and  update  the  density  fields. 

•  Solve  the  momentum  equations  explicitly  and  calculate  the  fields. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P,and 

•  Return  to  step  one  and  iterate  until  convergence _ 


The  MCBA  following  SIMPLEX  (MCBA-SIMPLEX):  Symbolic  Form 


Predictor: 


=  HP^[u®*  ]  -  pF 

Corrector: 

=  F  +F,p^'‘^*  =  p^^"  +p'*0 

:.vr  =HP^[u<''>**]-P*'”D^*VpP*  =  HP,[u<'‘>*  +^0 

■■\p('^y  =  cYF 

Condition: 


H  S 

*  [+Ap  .S]+  A^  |V*)°p(*>°u^*> 


Approximation: 

Neglect  App*'°p^*y^*^^s]  and  let 

uf'  =  HPp[u‘*^']-  r^^DfVpF  =  -r^*’°D^y^VpP' 


(81) 

(82) 

(83) 

(84) 


(85) 


(86) 

(87) 
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Assume  that  the  pressure  difference  local  to  the  velocity  is  representative  of  all  pressure  differences 
i.e.  thus: 


Approximate  Equation: 


X  |-  ]-  (r(*>  D'^>"^VP')s| 

k 

-  Pp _ yp  Pp  )  r>  L  A 


Q+App>y*>'t/^*^*]l- 


A  Global  MCBA-SIMPLEX  Iteration 


•  Solve  implicitly  for  using  the  old  pressure  and  density  fields. 

•  Calculate  the  fields. 

•  Solve  implicitly  for  the  fields. 

•  Solve  the  pressure  correction  equation  using  these  fields. 

•  Correct  P  and 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 

The  MCBA  following  SIMPLEM  (MCBA-SIMPLEM):  Symbolic  Form 

First  Predictor: 

No  predictor  stage.  Only  coefficients  of  the  momentum  equations  are  calculated. 

First  Corrector: 

+  =  P°  +P',p^^'’  =  +p^^0 

u<;>'  =  HP^[u^*^°]  -  VpP° 

Condition: 


rp-p^'- 


.£l  +  A,[r<‘V*’'U“’‘] 


=  -£■1  *  J-L  --  J  j. 

*  [+ A^  p^*’°  (hp[u^*^'])s]+  Ap p^*^'u®'.s]J 


Approximation: 


Neglect: 
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HP  ],  =5>  u^P'  = 

Approximate  Equation: 

Ap  VP')s]| 

jL  \_ijl  J 


=  _y  j IlIl — )...  a+^, 


Second  Predictor: 

u***  =HP,[u^*^**]-y^“DrV,P* 
Second  Corrector: 

No  corrector  stage. 

A  Global  MCBA-SIMPLEM  Iteration 


(95) 


(96) 


•  Calculate  the  fields  based  on  values  from  the  previous  iteration. 

•  Solve  the  pressure  correction  equation. 

•  Correct  P  and 

•  Calculate  new  and  fields. 

•  Solve  implicitly  for  using  the  new  fields. 

•  Solve  implicitly  for 

•  Solve  implicitly  the  energy  equations  and  update  the  density  fields. 

•  Return  to  the  first  step  and  iterate  until  convergence. _ 

The  Expanded  Form  of  the  Pressure-Correction  Equation 


In  simplified  form,  the  pressure-correction  equation  may  be  written  as: 

^  jik>C(k)p/ +  Ap  Ap  ®(r('')  D<'^^ 

K  LOl  ^  “* 

(97) 

Where,  depending  on  the  algorithm  used,  and  represent  values  from  the  previous  iteration  or 
from  a  previous  corrector  step.  The  second  term  on  the  left  hand  side  is  similar  to  a  convection  term 
and  may  be  discretized  using  any  convective  scheme  (the  upwind  scheme  is  employed  here).  The 

term  Ap  pkrp(k)^(kru(k)yp,^sj  is  discretized  following  the  same  procedure  used  in  discretizing 

the  diffusion  flux.  Substituting  the  various  terms  in  Eq.  (97)  by  their  equivalent  expressions  and 
neglecting  the  non-orthogonal  terms,  the  final  form  of  the  pressure-correction  equation  is  wntten  as: 

Arpp=SA:;p„',+Br 

nb 


-mT, 


5t 


•Q+App’y^v^)] 


(98) 
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where 

Aj-  =  S^“  +  (r“>-C<‘>)Jhu“,o|) 


Ar=iA;'+xi 

NB 

Br=-i| 


f  (r"'C**),Q> 

5t 


+ 


_  -(k)old  (k)old\ 

V^P  Pp  Tp  Pp  _l  V 


8t 

l(k)r 


i(r<"-crK 

nb 


(k) 


(99) 

(100) 


(101) 


imuXSrO 

— — - 

Sf  •Qf 

The  corrections  are  then  applied  to  the  velocity,  pressure,  and  density  fields  using  the  following 
equations: 


u«‘  =  u®“  -  r^'^^D'p'^VpP',  P*  =  P”  +  P',  =  p"''“  +  C®P'  (102) 

As  detailed  in  [59],  numerical  experiments  using  the  above  approach  to  simulate  air-water  flows  have 

shown  poor  conservation  of  the  lighter  fluid.  Without  going  into  details,  this  problem  can  be 

considerably  alleviated  by  normalizing  the  individual  continuity  equations,  and  hence  the  global  mass 

conservation  equation,  by  means  of  a  weighting  factor  such  as  a  reference  density  (which  is  fluid 

dependent).  This  approach  has  been  adopted  in  solving  all  problems  presented  in  this  work. 


Solving  the  energy  equations 

The  solutions  of  the  energy  equations  follow  that  of  the  general  multi-fluid  scalar  equation.  As  such, 
nothing  new  needs  to  be  added  in  that  regard  (though  in  many  cases  coupling  of  the  energy  equation 
with  the  momentum  and  continuity  equations  is  beneficial). 


Results  and  Discussion 


The  performance  of  the  various  multi-fluid  Mass  Conservation  Based  Algorithms  is  assessed  in  this 
section  by  presenting  solutions  to  several  one  and  two-dimensional  two-phase  problems.  Results  are 
presented  in  terms  of  the  CPU-time  needed  to  converge  the  solution  to  a  set  level  and  of  the 
convergence  history.  Moreover,  solutions  are  obtained  for  a  number  of  grids  in  order  to  assess  the 
performance  of  the  various  algorithms  with  increasing  grid  density.  For  the  two-dimensional 
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problems,  in  addition  to  the  CPU-time  needed  to  solve  a  problem  using  the  single  grid  method  (SG), 
the  CPU-time  needed  using  different  solution  strategies  is  also  displayed.  These  include  the 
prolongation  scheme  (PG)  in  which  the  solution  from  the  next  coarser  grid  is  used  to  provide  the 
initial  field,  and  the  full  multi-grid  method  (FMG).  Results  are  compared  against  available 
experimental  data  and/or  numerical/theoretical  values.  The  residual  of  a  variable  (|)  at  the  end  of  an 
outer  iteration  is  defined  as: 

RESf^  =  XAp(l)J^-  XA„b<t)lV-Bf  (103) 

c.v  all  p  neighbours 

For  global  mass  conservation,  the  imbalance  in  mass  is  defined  as: 

RES ,  =  XE  >  Q  -  Ap [r®p^'' (104) 

k  c.v. 

All  residuals  are  normalized  by  their  respective  inlet  fluxes.  Computations  are  terminated  when  the 
maximum  normalized  residual  of  all  variables,  drops  below  a  very  small  number  fig.  For  a  given 
problem,  the  same  value  of  £5  is  used  with  all  algorithms.  In  general,  it  is  found  that  requiring  the 
overall  mass  residuals  to  be  satisfied  to  within  Eg  is  a  very  stringent  requirement  and  the  last  to  be 
fulfilled.  This  is  why  these  residuals  are  the  ones  presented  here  and  used  to  compare  the  performance 
of  the  various  algorithms.  In  all  problems,  the  first  phase  represents  the  continuous  phase  (denoted  by 
a  superscript  (c)),  which  must  be  fluid,  and  the  second  phase  is  the  disperse  phase  (denoted  by  a 
superscript  (d)),  which  may  be  solid  or  fluid.  Unless  otherwise  specified  the  HR  SMART  scheme  is 
used  in  all  computations  reported  in  this  study.  For  a  given  problem,  all  results  are  generated  starting 
from  the  same  initial  guess.  Moreover,  it  should  be  stated  that  in  iterative  techniques,  different  initial 
guesses  might  require  different  computational  efforts. 

One-dimensional  two-phase  validation  problems 

Due  to  the  large  number  of  parameters  affecting  the  performance  of  the  various  algorithms  and  to 
allow  a  thorough  testing  of  these  algorithms,  eight  one-dimensional  two-phase  problems  are 
considered.  These  problems  can  be  broadly  classified  as:  ^  horizontal  particle  transport,  and  (ii) 
vertical  particle  transport.  Despite  its  geometric  simplicity,  the  one  dimensional  particle  transport 
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problem  represents  a  wide  range  of  physical  conditions.  The  effects  of  grid  refinement  on  accuracy 
and  convergence  are  studied  by  solving  the  problems  on  four  grid  systems  of  sizes  20, 40,  80,  and  160 
control  volumes  with  8s  assigned  the  value  of  10'^ 

Many  runs  were  performed  so  as  to  set  the  control  parameters  of  each  algorithm  near  optimum  values. 
The  task  was  performed  once  and  was  not  repeated  for  the  2-dimensional  problems.  The  CPU  times 
are  reported  in  the  form  of  charts  and  not  via  numerical  values.  This  allows  a  comparative  assessment 
of  their  performance.  Moreover,  all  CPU  times  are  normalized  by  the  time  needed  by  MCBA- 
SIMPLE  to  reach  the  set  residuals  on  the  coarsest  grid. 

Horizontal  particle  transport 

The  physical  situation  is  depicted  in  Fig.  2.  Depending  on  the  set  densities,  it  represents  either  the 
steady  flow  of  solid  particles  suspended  in  a  fi*ee  stream  of  air  or  the  steady  flow  of  air  bubbles  in  a 
stream  of  water.  The  slip  between  the  phases  determines  the  drag,  which  is  the  sole  driving  force  for 
the  particle-bubble/air-water  motion  (g=0).  In  the  suspension,  the  inter-particle/bubble  forces  are 
neglected.  Diffusion  in  both  phases  is  set  to  zero  and  the  inter-phase  drag  force  is  calculated  as: 


«  Tp 


8^. 


(105) 

(106) 

(107) 

The  drag  coefficient,  Q),  is  set  to  0.44.  Since  diffusion  is  neglected,  the  MCBA-SIMPLEST  and 
MCBA-PRIME  becomes  identical  and  reference  will  be  made  to  MCBA-SIMPLEST  only.  The  task  is 
to  calculate  the  particle/bubble-velocity  distribution  as  a  function  of  position.  If  the  flow  field  is 
extended  far  enough  (here  computations  are  performed  over  a  length  of  L=2m),  the  particle/bubble 
and  fluid  phases  are  expected  to  approach  an  equilibrium  velocity  given  by: 


^equilibrium  inlet  inlet  inlet  inlet 


(108) 


Dilute  gas-solid  flow 

The  steady  flow  of  dilute  particles  suspended  in  a  free  stream  of  air  is  studied  first.  At  inlet,  the  air 
and  particle  velocities  are  5  m/s  and  1  m/s,  respectively.  The  physical  properties  of  the  two  phases 
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are:  =  2000,  =  10  ^  Tp  =  1  mm,  =  10  ^  Due  to  the  dilute  concentration  of  the 

particles,  the  free  stream  velocity  is  more  or  less  unaffected  by  their  presence  and  the  equilibrium 
velocity  is  nearly  equal  to  the  inlet  free  stream  velocity.  Based  on  this  observation,  Morsi  and 
Alexander  [80]  obtained  the  following  analytical  solution  for  the  particle  velocity  u***'  as  a  function  of 
the  position  x  and  the  properties  of  the  two  phases: 


ytc) 

^  inlet 


V(<=) 

^  inlet  inlet 


(109) 


_ YmlsJ _ =  — £ _ .^x  +  Lnfv^'^  1- 

^L'' inlet  ^  S  r  ^  ^ L'' Met  '' inlet  J 

This  case  is  of  particiilar  importance  since  the  flow  situation  has  an  exact  solution.  As  shown  in  Fig. 
3(a)  the  predicted  particle  velocity  distribution  falls  on  top  of  the  analytical  solution  given  by  Eq. 
(109),  which  is  an  indication  of  the  accuracy  of  the  numerical  procedure.  The  convergence  history  of 
the  various  MCBA  and  over  the  four  grid  networks  used  are  displayed  in  Figs.  3(b)-3(h).  For  all 
algorithms,  the  required  number  of  iterations  increases  as  the  grid  size  is  increased,  with  PISO  (Fig. 
3(b))  requiring  the  minimum  and  SIMPLEST/PRIME  (Fig.  3(f))  the  maximum  number  of  iterations 
on  all  grids.  The  convergence  histories  of  SIMPLE,  SIMPLEC,  and  SIMPLEX  (Figs.  3(c),  3(d),  and 
3(g),  respectively)  are  very  similar  requiring  nearly  the  same  number  of  iterations  on  all  grids.  The 
convergence  path  of  SIMPLEM  (Fig.  3(e))  is  not  as  smooth  as  that  of  SIMPLE  due  to  the  fact  that  at 
the  end  of  an  outer  iteration,  the  velocity  field  is  momentum  satisfying  rather  than  continuity 
satisfying.  The  convergence  paths  of  the  various  algorithms  over  a  grid  of  size  80  C.V.  are  compared 
in  Fig.  3(h)  and  the  above  observations  can  easily  be  inferred  from  the  figure.  The  normalized  CPU 
efforts  required  by  the  various  algorithms  over  all  grids  are  depicted  in  Fig.  4.  As  expected,  the  chart 
clearly  shows  that  the  CPU  time  increases  with  increasing  grid  density.  Moreover,  it  is  hard  to  see  any 
noticeable  difference  in  the  CPU  times  for  SIMPLE,  SIMPLEC,  SIMPLEX,  and  PISO.  The 
SIMPLEM  algorithm  requires  slightly  higher  computational  effort  as  compared  to  SIMPLE.  The 
worst  performance  is  for  SIMPLEST  which  degenerates  to  PRIME  in  the  absence  of  diffusion  and 
results  in  a  fully  explicit  solution  scheme. 


Dense  gas-solid  flow 
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The  only  difference  between  this  case  and  the  previous  one  is  in  the  concentration  of  particles,  which 
is  set  to  =  10”^ .  Despite  the  low  value  of  the  inlet  disperse  phase  volume  fraction,  this  problem  is 
dense  in  the  sense  that  the  ratio  of  disperse  phase  and  continuous  phase  mass  loadings  is  large 
r(d)p(d)/r<^>p<"^  =  20.  Thus  the  disperse  phase  carries  most  of  the  inertia  of  the  mixture.  The 
equilibrium  velocity  in  this  case,  as  obtained  from  Eq.  (108)  is  4.96  m/s  as  compared  to  4.99996  m/s 
in  the  previous  case.  Due  to  this  slight  difference  between  the  inlet  air  velocity  and  the  final 
equilibrium  velocity,  the  free  stream  velocity  may  be  assumed  to  be  nearly  constant  and  the  variation 
in  particle  velocity  can  be  obtained  from  Eq.  (109)  (i.e.  the  variation  in  particle  velocity  is  exactly  the 
same  as  in  the  previous  case).  The  predicted  air  and  particle  velocity  distributions  are  displayed  in 
Fig.  5(a).  The  numerical  and  anal34ical  particle  velocity  profiles  are  indistinguishable  and  fall  on  top 
of  each  others  (denoted  solid  in  the  Figure).  Moreover,  the  slight  decrease  in  the  air  velocity  can  be 
easily  depicted.  The  convergence  paths  for  all  algorithms  and  over  all  grid  systems  used  are  displayed 
in  Figs.  5(b)-5(h).  In  general,  a  larger  number  of  iterations  are  required  to  reach  the  desired  level  of 
convergence  on  a  given  grid  as  compared  to  the  dilute  case  due  to  the  increased  importance  of  the 
inter-phase  term.  The  general  convergence  trend  is  similar  to  that  of  the  dilute  problem  with  PISO 
requiring  the  minimum  and  SIMPLEST  the  maximum  number  of  iterations.  The  SIMPLEX  algorithm 
(Fig.  5(g))  is  seen  to  require  a  slightly  lower  number  of  iterations  on  the  finest  grid  as  compared  to 
SIMPLE  (Fig.  5(c)),  SIMPLEC  (Fig.  5(d)),  and  SIMPLEM  (Fig.  5(e)).  The  smoothest  convergence 
paths  are  for  SIMPLEX  and  SIMPLEC.  As  depicted  in  Figs.  5(f)  and  5(h),  the  performance  of 
SIMPLEST/PRIME  is  poor  as  compared  to  other  algorithms  for  the  same  reasons  stated  above.  The 
normalized  CPU  times  required  by  the  various  algorithms  on  all  grids  are  presented  in  Fig.  6.  The 
same  conclusions  drawn  above  are  applicable  here  with  PISO  requiring  the  lowest  computational 
effort  (10%  less  than  SIMPLE  on  the  finest  mesh).  The  performance  of  SIMPLE,  SIMPLEC,  and 
SIMPLEX  is  more  or  less  identical  while  that  of  SIMPLEM  is  of  lower  quality  necessitating 
increasingly  higher  computational  effort  with  denser  meshes,  and  requiring  about  40%  more  effort 
than  SIMPLE  on  the  fine  grids  (80  and  160  C.V.).  The  computational  effort  needed  by 
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SIMPLEST/PRIME  is  however  the  most  extensive  and  is  nearly  623%  the  one  needed  by  SIMPLE  on 
the  finest  mesh. 

Dilute  bubbly  flow 

For  the  same  configuration  displayed  in  Fig.  2,  the  continuous  phase  is  considered  to  be  water  and  the 
disperse  phase  to  be  air.  The  resulting  flow  is  denoted  in  the  literature  by  bubbly  flow.  With  the 
exception  of  =  10”^  and  at  inlet  =0.1,  other  physical  properties  and  inlet  conditions  are 

the  same  as  those  considered  earlier.  This  is  a  strongly  coupled  problem  and  represents  a  good  test  for 
the  numerical  procedure  and  performance  of  the  algorithms.  The  correct  physical  solution  is  that  the 
bubble  and  continuous  phase  velocities  both  reach  the  equilibrium  velocity  of  4.6  m/s  (Eq.  (108))  in  a 
distance  too  small  to  be  correctly  resolved  by  any  of  the  grid  netwroks  used.  Results  for  this  case  are 
presented  in  Figs.  7  and  8.  Axial  velocity  distribution  for  both  water  and  air  are  displayed  in  Fig.  7(a). 
Both  phases  reach  the  equilibrium  velocity  of  4.6  m/s  over  a  very  short  distance  from  the  inlet  section 
and  remain  constant  afterward,  as  expected.  The  relative  convergence  characteristics  of  the  various 
algorithms  remain  the  same.  However,  all  algorithms  require  larger  number  of  iterations  as  compared 
to  the  dilute  gas  solid  flow  case  due  to  the  stronger  coupling  between  the  phases.  Consistently,  the 
PISO  (Fig.  7(b))  and  SIMPLEST/PRIME  (Fig.  7(f))  algorithms  need  the  lowest  and  highest  number 
of  iterations,  repsectively.  As  in  the  previous  two  cases,  the  convergence  attributes  of  SIMPLE  (Fig. 
7(c)),  SIMPLEC  (Fig.  7(d)),  and  SIMPLEX  (Fig.  7(g))  are  very  similar  and  those  of  SIMPLEM  (Fig. 
7(e))  are  close  to  them.  The  large  difference  in  performance  between  SIMPLEST/PRIME  and  the 
remaining  algorithms  is  clearly  demonstarted  in  Fig.  7(h).  The  normalized  CPU  time  of 
SIMPLEST/PRIME  for  this  problem  (Fig.  8)  is  lower  than  in  the  previous  two  problems  due  to  a 
higher  rate  of  increase  in  the  time  needed  by  other  algorithms  (the  real  computational  time  of  all 
algorithms  has  increased).  The  relative  performance  of  the  various  algorithms  is  nearly  as  described 
earlier  with  the  time  required  by  of  PISO,  SIMPLE,  SIMPLEC,  and  SIMPLEX  being  on  average  the 
same.  The  SIMPLEST/PRIME  algorithm  however,  requires  nearly  three  folds  the  time  needed  by 
SIMPLE,  which  represents  a  noticeable  improvement. 


Dense  bubbly  flow 
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The  only  difference  between  this  case  and  the  previous  one  is  in  the  concentration  of  bubbles,  which 
is  set  to  =0.5  .  With  such  high  value  of  void  fraction,  bubble  coalescence  may  occur.  However, 
this  is  not  accounted  for  here.  The  analytical  solution  is  the  same  as  in  the  previous  case  with  the 
equilibrium  velocity,  as  computed  from  Eq.  (108),  being  3  m/s.  As  depicted  in  Fig.  9(a),  the 
equilibrium  velocity  obtained  numerically  is  exact.  Moreover,  the  performance  of  the  various 
algorithms  (Figs.  9(b)-9(h))  and  their  normalized  CPU  times  (Fig.  10)  vary  relatively  in  a  manner 
similar  to  what  was  previously  discussed  and  it  is  deemed  redundant  to  be  repeated. 


Vertical  particle  transport 

Here,  the  flow  is  in  the  vertical  direction  (Fig.  1)  and  the  gravitational  acceleration  is  assigned  the 
constant  value  of  g=10  m/s^  and  the  flow  field  is  extended  a  length  of  L=20m.  For  this  situation,  the 
velocities  of  the  two  phases  do  not  equilibrate  to  a  constant.  Rather,  the  disperse  phase  equilibrates  to 
a  finite  settling  velocity  relative  to  the  continuous  phase,  at  which  the  gravitational  force  balances  the 
drag  force  [81].  As  for  the  horizontal  transport  problems,  the  inter-particle/bubble  forces  are 
neglected.  However,  unlike  the  previous  situation,  diffusion  in  the  continuous  phase  is  retained. 
Moreover,  the  inter-phase  drag  force  is  calculated  using  Eqs.  (105)-(106)  and  the  drag  coefficient,  Cd, 
is  considered  to  be  particle  Reynolds  number  dependent  and  calculated  as: 


Cd  =  — +0.44 

Rs 

where 


(110) 


2r  V 

^lip 


(111) 


Since  diffusion  in  the  continuous  phase  is  not  neglected,  the  performances  of  MCBA-SIMPLEST  and 
MCBA-PRIME  are  expected  to  be  different. 


Dilute  gas-solid  flow 

The  material  properties  and  boundary  conditions  considered  for  this  case  are  given  by: 

p(‘‘>/p('>  =  1000,  1^=1  mm 

Vl;:,i=100  m/s,  VL‘;>,  =  10  m/s, 


(112) 

(113) 
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The  large  velocity  boundary  conditions  are  used  to  ensure  that  the  solid  phase  does  not  exit  the  inlet. 
The  predicted  air  and  particle  velocity  distributions  depicted  in  Fig.  11(a)  are  in  excellent  agreement 
with  similar  predictions  reported  in  [81].  As  shown  in  Figs.  11(b)- 11(h),  the  decrease  in  the  mass 
residuals  is  highly  non-monotonic  showing  a  cyclic  decaying  behavior.  The  PISO  algorithm  (Fig. 
11(b))  seems  to  be  the  least  affected  and  its  convergence  path  shows  little  cycling  as  compared  to 
other  algorithms  and  requiring  the  least  number  of  iterations  over  all  grids.  It  should  be  mentioned 
here  that  decreasing  the  under-relaxation  factors  could  have  smoothed  this  cycling  behavior  out.  This, 
however,  would  have  been  accomplished  at  the  expense  of  increasing  the  computational  time.  Even 
though  the  convergence  is  not  monotonic,  the  number  of  iterations  needed  to  converge  the  solution  to 
the  desired  level  is  very  close  to  that  needed  in  the  similar  horizontal  transport  case.  The  convergence 
of  the  SIMPLEC  algorithm  over  the  dense  grid  system  (Fig.  2(d))  shows  cycles  of  large  amplitudes 
due  to  the  fact  that  pressure  is  not  under-relaxed.  The  performance  of  SIMPLEST  (Fig.  2(f))  and 
PRIME  (Fig.  2(g))  is  veiy  close,  with  SIMPLEST  requiring  slightly  lower  number  of  iterations  due  to 
the  small  implicitness  introduced  by  the  diffusion  of  the  continuous  phase.  However,  both  require  on 
the  finest  mesh  almost  1300%  the  number  of  iterations  needed  by  SIMPLE.  The  number  of  outer 
iterations  needed  by  SIMPLEX,  SIMPLEC,  and  SIMPLEM  is  very  close  to  that  of  SIMPLE,  with 
SIMPLEM  requiring  the  highest  number  of  iterations.  This  behavior  is  further  revealed  by  the 
normalized  CPU  times  needed  by  the  various  algorithms  that  are  displayed  in  Fig.  12.  The  efficiency 
of  SIMPLEST  is  slightly  better  than  PRIME,  both  however  are  about  five  times  more  expensive  than 
all  other  algorithms  whose  performance  is  very  comparable. 

Dense  gas-solid  flows 

The  material  properties  and  boundary  conditions  are  similar  to  the  previous  case  with  the  exception  of 
the  particles’  volume  fraction,  which  is  set  to  =  10~^.  Predicted  air  and  particle  velocity  profiles 
are  displayed  in  Fig.  13(a)  while  mass  residuals  are  presented  in  Figs.  13(b)-13(h).  The  convergence 
paths  are  smoother  in  comparison  with  the  dilute  case  due  to  the  higher  volume  fraction  of  the 
disperse  phase,  which  makes  the  computations  less  sensitive  to  the  intermediate  level  of  convergence. 
However,  higher  number  of  iterations  is  needed  in  comparison  with  the  dilute  case  due  to  the  higher 
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mass-loading  ratio.  Besides  that,  the  convergence  behavior  is  similar  to  the  previous  cases  with 
SIMPLEST  (Fig.  13(f))  and  PRIME  (Fig.  13(g))  requiring  the  highest  number  of  iterations 
(SIMPLEST  needs  slightly  less  iterations  for  the  reason  stated  above).  The  number  of  iterations 
needed  by  PISO,  SIMPLEC,  and  SIMPLEX  (compare  Figs.  13(b),  13(d),  and  13(h))  is  very  close. 
The  SIMPLEM  algorithm  (Fig.  13(e))  takes  a  slightly  higher  number  of  iterations  than  SIMPLE  (Fig. 
13(c)).  For  this  problem  SIMPLEX  requires  the  least  computational  time  followed  by  SIMPLEC  and 
SIMPLE  (Fig.  14).  Again  the  computational  effort  needed  by  SIMPLEST  and  PRIME  is  much  higher 
(about  700%  the  time  needed  by  SIMPLEX  on  the  finest  grid). 

Dilute  bubbly  flows 

In  this  problem,  the  continuous  phase  is  water  and  the  disperse  phase  is  air.  With  the  exception  of 
set  to  10"^  and  to  0.1  at  inlet,  other  physical  properties  and  inlet  conditions  are  the  same 

as  those  considered  earlier.  This  is  a  very  difficult  problem  to  get  convergence  to  unless  the  proper 
under-relaxation  is  used.  It  was  possible  to  get  feasible  solutions  (i.e  with  reasonable  computational 
time)  when  under-relaxing  by  inertia  (i.e.  through  the  use  of  false  time  steps).  For  the  results 
presented  in  Figs.  15  and  16,  a  time  step  (At)  of  value  10""*  s  is  used  for  the  velocity  field  of  the 
dispersed  gas  phase,  At=l  s  for  the  volume  fractions,  and  At=0.01  s  for  the  velocity  field  of  the  liquid 
phase  (this  last  value  is  employed  with  all  algorithms  except  PISO  and  PRIME  for  which  a  value  of 
0.1  s  is  used).  Also,  to  accelerate  convergence,  it  is  found  advantageous  not  to  under-relax  the 
pressure  field.  As  expected,  the  correct  velocity  fields  are  predicted  (Fig.  15(a)).  All  algorithms 
require  higher  number  of  iterations  as  compared  to  the  previous  problems.  Moreover,  convergence 
flattens  after  a  certain  level  because  of  the  large  under-relaxation  used  for  the  bubble  phase 
momentum  equations.  Using  an  automatic  time  stepping  algorithm  that  increases  this  time  step  as  the 
residuals  are  decreased  could  resolve  this  issue.  The  performance  of  PISO  (Fig.  15(b))  is  totally 
unexpected  requiring  higher  number  of  iterations  than  SIMPLE  (Fig.  15(c)),  SIMPLEC  (Fig.  15(d)), 
SIMPLEM  (Fig.  15(e)),  and  SIMPLEX  (Fig.  15(h)).  The  decrease  in  the  rate  of  convergence  as  the 
residuals  are  decreased,  is  very  clear  with  SIMPLE,  SIMPLEC,  and  SIMPLEX  which  show  a  sudden 
change  in  their  convergence  slopes.  The  convergence  histories  of  these  algorithms  are  nearly 
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identical.  The  explicitness  introduced  in  PISO  (Fig.  15(b)),  SIMPLEST  (Fig.  15(f)),  and  PRIME  (Fig. 
15(g))  causes  an  increase  in  the  number  of  iterations  and  a  highly  non-monotonic  convergence 
behavior  and  seems  to  be  undesirable  in  laminar  bubbly  flows.  The  rate  of  convergence  of  the 
SIMPLEM  algorithm  (Fig.  15(e))  is  nearly  constant  and  does  not  show  any  change  in  slope  as  the 
computations  progress.  This  is  due  to  the  larger  number  of  inner  iterations  needed  on  the  pressure 
correction  equation  to  promote  convergence  and  is  associated  with  higher  computational  cost  as 
depicted  in  Fig.  16.  The  noticeable  change  in  the  normalized  time  chart  is  the  performance  of 
SIMPLEST  which  shows  a  good  improvement  and  that  of  PISO,  which  has  deteriorated.  The  CPU 
times  needed  by  SIMPLE,  SIMPLEC,  and  SIMPLEX  are  consistently  very  close.  For  a  given  grid 
however,  the  largest  CPU  time  consumed  is  only  double  the  lowest  CPU  needed. 

Dense  bubbly  flows 

With  the  exception  of  setting  to  0.5,  the  physical  situation,  material  properties,  and  boundary 
conditions  are  the  same  as  in  the  previous  problem.  Results  for  the  problem  are  presented  in  Figs.  17 
and  18.  The  predicted  liquid  and  gas  velocity  distributions,  which  are  in  excellent  accord  with 
published  data,  are  depicted  in  Fig.  17(a).  The  trend  of  convergence  (Figs.  17(b)-17(h))  is  very  similar 
to  the  dilute  case  with  the  following  differences.  With  PISO  and  PRIME,  in  order  to  drive  residuals  to 
the  desired  level  of  convergence,  the  SMART  scheme  had  to  be  blended  with  the  UPWIND  scheme. 
The  percentage  increased  from  15%  on  the  coarsest  grid  to  50%  on  the  finest  grid.  Moreover,  mass 
residuals  reached  the  desired  convergence  level  long  before  the  momentum  residuals.  Thus,  the 
numbers  of  iterations  needed  by  PISO  and  PRIME  are  much  higher  than  the  ones  presented  in  Figs. 
17(b)  and  17(g)  and  are  reflected  by  the  CPU  times  displayed  in  Fig.  18.  In  addition,  to  stabilize 
SIMPLEST  on  the  finest  grid,  the  SMART  scheme  had  to  be  blended  with  20%  of  the  upwind  value. 
The  difficulties  faced  by  PISO,  PRIME,  and  SIMPLEST  are  attributed  to  the  additional  explicitness 
introduced  as  a  result  of  using  a  HR  scheme  implemented  via  a  deferred  correction  strategy  [82].  This 
deferred  correction  technique  has  lower  influence  on  the  performance  of  the  remaining  algorithms  due 
to  their  higher  implicitness.  As  shown  in  Fig.  18,  even  with  the  blending  strategy,  PISO  and  PRIME 
still  need  a  high  computational  time  to  decrease  the  residuals  to  the  desired  level.  The  SIMPLEST 
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algorithm  requires  the  lowest  computational  effort  among  the  three  algorithms  on  the  finest  grid  when 
blended  with  the  upwind  scheme.  The  CPU  times  consumed  by  SIMPLE,  SIMPLEC,  and  SIMPLEX 
are  consistently  very  close  to  each  other. 

By  comparing  the  behavior  of  the  various  algorithms  in  all  problems,  it  is  clear  that  the  performance 
of  SIMPLE,  SIMPLEC,  and  SIMPLEX  is  consistent  and  require,  on  average,  the  least  computational 
effort.  Even  though  in  some  cases  PISO  consumes  less  CPU  time,  its  performance  for  upward  bubbly 
flows  is  not  satisfactory.  The  performance  of  SIMPLEM  is  consistent,  but  demands  more  CPU  time 
than  SIMPLE,  SIMPLEC,  and  SIMPLEX.  The  SIMPLEST  algorithm  is  recommended  for  diffusion- 
dominated  flows  (e.g.  turbulent  flows).  Nevertheless,  other  algorithms  may  perform,  in  such 
situations,  equally  well  if  not  better.  The  PRIME  algorithm  is  the  most  expensive  to  use  on  all  grids 
and  for  all  physical  situations  presented  here.  Most  importantly  however,  is  the  fact  that  all  these 
algorithms  can  be  used  to  predict  multi-phase  (in  this  case  two-phase)  flows.  The  next  step  is  to  test 
these  algorithms  in  multi-dimensional  incompressible,  compressible,  laminar,  and  turbulent  flows. 

Two-dimensional  two-phase  validation  problems 

In  this  section,  four  two-dimensional  two-phase  problems  are  solved.  The  first  two  problems  deal  with 
incompressible  turbulent  flows  while  the  last  two  problems  are  concerned  with  compressible  flows. 
As  mentioned  earlier,  solutions  will  be  obtained  using  the  prolongation  grid  scheme  (PG)  and  the  full 
multi-grid  method  (FMG)  in  addition  to  the  usual  single  grid  (SG)  approach.  This  will  enable 
assessing  the  performance  of  the  non-linear  multi-grid  with  multi-phase  flows. 

Turbulent  upward  bubbly  flow  in  a  pipe 

The  problem  considered  involves  the  prediction  of  radial  phase  distribution  in  turbulent  upward  air- 
water  flow  in  a  pipe.  Many  experimental  and  numerical  studies  addressing  this  problem  have 
appeared  in  the  literature  [83-91].  Most  of  these  studies  have  indicated  that  the  lateral  forces  that  most 
strongly  affect  the  void  distribution  are  the  turbulent  stresses  and  the  lateral  lift  force.  As  such,  in 
addition  to  the  usual  drag  force,  the  lift  force  is  considered  as  part  of  the  interfacial  force  terms  in  the 
momentum  equations.  In  the  present  work,  the  interfacial  drag  forces  per  unit  volume  are  given  by: 
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Cl  =  0,3(1-  2.78(0.2,  (119) 

where  (a,  b)  denotes  the  minimum  of  a  and  b  and  Cla  is  an  empirical  constant. 

Besides  the  drag  and  lift  interfacial  forces,  the  effect  of  bubbles  on  the  turbulent  field  is  very 

important,  since  the  distribution  of  bubbles  affects  the  turbulence  field  in  the  liquid  phase  and  at  the 

same  time  the  liquid  phase’s  turbulence  is  influenced  by  the  bubbles.  In  this  work,  turbulence  is 

assumed  to  be  a  property  of  the  continuous  liquid  phase  (c)  and  the  turbulent  kinematic  viscosity  of 
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the  dispersed  air  phase  (d)  is  assumed  to  be  a  function  of  that  of  the  continuous  phase.  The  turbulent 
viscosity  of  the  continuous  phase  is  computed  by  solving  the  following  modified  transport  equations 
for  the  turbulent  kinetic  energy  k  and  its  dissipation  rate  e  that  take  into  account  the  interaction 
between  the  phases; 
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where  is  the  well  known  volumetric  production  rate  of  1^^^  by  shear  forces,  the  turbulent 
Schmidt  number  for  volume  fractions,  and  Pb  is  the  production  rate  of  by  drag  due  to  the  motion 
of  the  bubbles  through  the  liquid  and  is  given  by: 

P  ,  0^15C,C,,'‘VV%, 

■p 

In  Eq.  (122)  Cb  is  an  empirical  constant  representing  the  fraction  of  turbulence  induced  by  bubbles 
that  goes  into  large-scale  turbulence  of  the  liquid  phase.  Moreover,  as  suggested  in  [91],  the  flux 
representing  the  interaction  between  the  fluctuating  velocity  and  volume  fraction  is  modeled  via  a 

gradient  diffusion  approximation  and  added  as  a  source  term  in  the  continuity  (v.^p^^^^D^’^Vr^^^^^and 
momentum  equations  with  the  diffusion  coefficient  D  given  by: 


The  turbulent  kinematic  eddy  viscosity  of  the  dispersed  and  continuous  phases  are  related  through: 

■^(c) 

(124) 

where  CJf  is  the  turbulent  Schmidt  number  for  the  interaction  between  the  two  phases.  The  above 


described  turbulence  model  is  not  exactly  the  one  described  in  [91]  but  rather  a  modified  version  of 
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that  model  in  which  the  turbulent  kinematic  viscosities  of  both  phases  are  allowed  to  be  different  in 
contrast  to  what  is  done  in  [91].  This  is  accomplished  through  the  introduction  of  the  Of  parameter. 
As  such,  different  diffusion  coefficients  are  used  for  the  different  phases. 

To  check  the  validity  and  correctness  of  the  above  described  treatment  two  different  situations,  for 
which  data  are  available,  are  considered.  In  the  first,  the  data  reported  by  Seriwaza  et  al  [83]  is 
reproduced.  The  second  situation  however,  is  concerned  with  the  data  of  Lahey  et  al  [94].  For  the 
Seriwaza  et  al  data  [83],  the  Reynolds  number  based  on  superficial  liquid  velocity  and  pipe  diameter 
is  8x10^,  the  inlet  superficial  gas  and  liquid  velocities  are  0.077  and  1.36  m/s,  respectively,  and  the 
inlet  void  fraction  is  5.36x10'^  with  no  slip  between  the  incoming  phases.  Moreover,  the  bubble 
diameter  is  taken  as  3  mm  [91],  while  the  fluid  properties  are  taken  as  p^‘'^=1000  Kg/m^,  p^'*^=1.23 

Kg/m^  and  10"^  m^/s.  On  the  other  hand,  Lahey  et  al’s  data  [94]  is  for  a  Reynolds  number  based 
on  superficial  liquid  velocity  and  pipe  diameter  of  5x10^,  the  inlet  superficial  gas  and  liquid  velocities 
are  0.1  and  1.08  m/s,  respectively.  Both  problems  are  solved  using  the  same  values  for  all  constants  in 
the  model  with  Cia=0.075,  apO.5,  ar=0.7,  and  Q-O.OS.  Predicted  radial  profiles  of  the  vertical  liquid 
velocity  and  void  fraction  presented  in  Figs.  19(a)  and  19(b)  using  a  grid  of  size  96x32  control 
volumes  concur  very  well  with  measurements  and  compare  favorably  with  numerical  profiles  reported 
by  Boisson  and  Malin  [91]  (Fig.  19(a))  and  PHOENICS  [96]  (Fig.  19(b)). 

Having  established  the  credibility  of  the  physical  model  and  numerical  procedure,  the  next  task  is  to 
assess  the  merits  of  the  various  algorithms  for  such  flows.  For  that  purpose  the  Seriwaza  et  al  [83] 
configuration  is  chosen  and  the  axi-symmetric  two-dimensional  calculations  are  performed  on  4 
different  grids  of  sizes  12x4,  24x8,  48x16,  and  96x32  control  volumes.  On  each  grid,  solutions  are 
generated  using  the  SG,  PG,  and  FMG  strategies  for  all  algorithms.  Results  are  displayed  in  the  form 
of  $)  total  mass  residuals  summed  over  both  phases  as  a  function  of  outer  iterations,  and  (ii) 
normalized  CPU  time  (Table  1)  needed  for  the  maximum  normalized  residual  of  all  variables  and  for 
all  phases  to  drop  below  Es=10'^.  In  order  not  to  overload  plots,  only  residuals  over  the  densest  grid 
using  the  SG,  PG,  and  FMG  methodologies  are  presented  in  Fig.  20.  Moreover,  it  should  be  clarified 
that  in  the  case  of  PISO  and  PRIME  the  number  of  iterations  needed  to  satisfy  the  convergence 
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criteria  is  much  higher  than  the  one  presented  (this  will  be  reflected  by  their  respective  normalized 
CPU  times  displayed  in  Table  1)  due  to  the  slower  convergence  rate  of  the  axial  liquid  momentum 
and  turbulence  kinetic  energy  equations,  i.e.  the  convergence  rate  of  these  equations  is  much  slower 
than  the  overall  mass  conservation  equation.  Since  the  overall  mass  residxxals  are  presented,  only  the 
iterations  needed  to  drive  these  residuals  to  the  required  level  (10‘^)  are  displayed,  even  though  at  the 
state  of  convergence  these  residuals  are  much  lower  (of  the  order  of  10'^). 

As  can  be  seen  from  Fig.  20,  it  is  possible  to  converge  the  solution  to  the  desired  level  with  all 
algorithms.  Moreover,  the  performance  of  SIMPLE  (Fig.  20(b)),  SIMPLEC  (Fig.  20(c)),  SIMPLEST 
(Fig.  20(e)),  and  SIMPLEX  (Fig.  20(g))  is  very  close.  This  is  also  demonstrated  in  Fig.  20(h)  where 
the  convergence  histories  using  the  FMG  method  are  plotted  for  all  algorithms.  As  shown,  the  plots 
for  these  algorithms  are  hard  to  distinguish.  For  all  algorithms,  the  SG  requires  the  largest  number  of 
iterations  and  the  FMG  the  lowest  number  of  iterations.  The  FMG  method  greatly  reduces  the  number 
of  iterations  whereas  reduction  with  the  PG  method  is  not  significant.  The  normalized  CPU  times 
displayed  in  Table  1  further  reveals  this.  As  shown  in  the  plots,  the  rate  of  convergence  of  the  SG  and 
PG  methods  diminishes  as  the  iterations  proceed.  This  is  in  contrast  with  the  convergence  behavior  of 
the  FMG  method,  which  seems  to  maintain,  more  or  less,  the  same  convergence  rate  (on  the  finest 
level)  as  the  iterations  progress. 

In  Table  1,  the  normalized  CPU-tknes  required  by  the  different  algorithms  to  solve  the  problem  on  the 
various  grids  and  with  the  various  methodologies  are  presented.  In  addition,  the  ratio  of  the  time 
needed  by  the  SG  method  to  the  one  needed  by  either  the  PG  or  the  FMG  method  is  displayed.  This 
allows  for  a  direct  quantitative  assessment  of  their  acceleration  rate.  As  depicted,  the  CPU  effort 
increases  for  all  methods  and  algorithms  with  increasing  the  grid  size.  The  CPU  times  of  SIMPLE, 
SIMPLEC,  SIMPLEX,  and  SIMPLEST  are  very  on  all  grids  and  for  all  methods  with  no  clear 
superiority  to  any  algorithm  over  the  others.  PISO  requires  the  least  computational  time  on  the 
coarsest  grid  used  (12x4  C.V.)  however,  as  the  grid  size  increases  the  picture  is  reversed  and  its 
performance  degrades  and  become  closer  to  that  of  PRIME  and  SIMPLEM.  The  use  of  the  PG 
method  on  the  relatively  coarse  grids  increases  the  computational  cost.  The  benefit  of  the  PG  method 
is  realized  on  the  relatively  dense  grids.  On  the  densest  grid,  the  savings  in  computational  effort  vary 
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from  11%  for  PRIME  to  28%  for  SIMPLEM  with  an  average  acceleration  rate  for  all  algorithms  of 
22%,  The  use  of  the  FMG  method  is  highly  beneficial  and  its  virtues  generally  appear  on  all  grids 
with  the  rate  of  acceleration  increasing  with  increasing  grid  density.  On  the  finest  grid,  the  rate  of 
acceleration  varies  from  144%  for  SIMPLEM  to  555%  with  PRIME  with  an  average  acceleration  rate 
for  all  algorithms  of  444%.  This  average  rate  of  acceleration  is  nearly  20  times  the  one  obtained  with 
the  PG  method.  This  indicates  that  the  significant  non-linearity  in  multi-fluid  flows  is  properly 
resolved  using  the  non-linear  multi  grid  method.  Moreover,  turbulence  could  be  the  reason  behind  the 
poor  performance  of  the  PG  method. 


Turbulent  air-particle  flow  in  a  vertical  pipe 

The  problem  consists  of  predicting  the  upward  flow  of  a  dilute  gas-solid  mixture  in  a  vertical  pipe 
using  the  various  algorithms  and  calculation  procedures.  For  these  simulations,  the  axi-symmetric 
form  of  the  gas  and  particulate  transport  equations  are  employed.  As  reported  in  several  studies  [97- 
99],  the  effects  of  interfacial  virtual  mass  and  lift  forces  are  small  and  may  be  neglected.  The 
controlling  interfacial  force  is  drag  and  the  closure  used  to  compute  it  in  the  present  study  is  due  to 
Harlow  and  Amsden  [100]  in  which  the  viscous  part  is  neglected.  Denoting  the  continuous  and 
dispersed  phases  by  superscripts  (c)  and  (d),  respectively,  the  drag  in  the  x-  and  y-momentum 
equations  are  given  by: 


where  fp  represents  the  particle’s  radius,  Cd  the  drag  coefficient  computed  from: 


(125) 

(126) 


for  RCp  <  1 

for  1  <  RCp  <  1000 
for  RCp  >  1000 


(127) 


and  RCp  the  Reynolds  number  based  on  the  particle  size  as  defined  in  Eq.(l  11). 
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As  before,  turbulence  is  assumed  to  be  a  property  of  the  continuous  gas  phase  (c)  and  the  turbulent 
kinematic  viscosity  of  the  dispersed  particle  phase  (d)  is  assumed  to  be  a  function  of  that  of  the 
continuous  phase.  Turbulence  modulation  due  to  the  presence  of  particles  is  predicted  using  a  two- 
phase  k-8  model.  Several  extensions  of  the  k-8  model  for  carrier-phase  turbulence  modulation  have 
been  proposed  in  the  literature  [70-76]  and  the  one  suggested  by  Chen  and  Wood  [72],  which 
introduces  additional  source  terms  into  the  turbulence  transport  equations,  is  adopted  here.  These 
source  terms  are  always  negative  and  act  to  reduce  k  and  8.  However,  depending  on  the  relative  extent 
of  reductions  in  k  and  8,  the  turbulent  viscosity  may  be  either  reduced  or  increased  by  the  presence  of 
particles.  Thus,  the  turbulent  viscosity  is  computed  by  solving  the  turbulence  transport  equations  (Eqs. 
(5)  and  (6))  for  the  continuous  phase  with  4^^  and  evaluated  using  the  following  relations 
suggested  by  Chen  and  Wood  [72]: 


(i.e.  the  interaction  terms  included  for  bubbly  flows  are  neglected  here),  respectively.  Figure  21  shows 


j  /cJl>m/c-jUU.3cu  J  y  uiutuKr  ivic-nuju  vvun  f  mi  iYH4,iit,-\Jt  lu  unuti^  jut  iviuiii-i  ' lUiu  ^  - luvv  ui  nn  Liyc-t^uo 


-TU 


the  fully  developed  gas  and  particles  mean  axial  velocity  profiles  generated  using  a  grid  of  size  96x40 
C.V.  It  is  evident  that  there  is  generally  a  very  good  agreement  between  the  predicted  and 
experimental  data  with  the  gas  velocity  being  slightly  over  predicted  and  the  particles  velocity  slightly 
under  predicted.  Moreover,  close  to  the  wall,  the  model  predictions  indicate  that  the  particles  have 
higher  velocities  than  the  gas,  which  is  in  accord  with  the  experimental  results  of  Tsuji  et  al  [97]. 

After  establishing  the  credibility  of  the  physical  model  and  numerical  procedure,  the  problem  is 
solved  over  four  different  grids  of  sizes  12x5,  24,xl0,  48x20,  and  96x40  control  volumes  using  the 
SIMPLE,  SIMPLEC,  SIMPLEX,  and  SIMPLEST  multi-fluid  algorithms  and  the  SG,  PG,  and  FMG 
solution  methods.  In  solving  he  problem,  it  is  noticed  that  the  initial  guess  greatly  affects  the 
convergence  history  and  time  required  to  converge  the  solution  to  the  desired  level.  Except  when 
solving  on  the  finest  mesh  with  the  SG  method,  the  initial  guess  used  for  the  velocity  field  is 
m/s.  The  use  of  this  initial  value  with  the  SG  method  on  the  finest  mesh  increased  the  CPU 
effort  needed  by  500%  over  the  one  needed  when  starting  with  an  initial  field  of  u^^W^‘^^=15.6  m/s.  To 
reduce  cost,  the  latter  initial  guess  is  used.  For  this  reason  the  mass  residuals  start  from  somehow  a 
lower  value  than  expected.  As  in  the  previous  problem,  results  are  displayed  in  the  form  of  $)  total 
mass  residuals  summed  over  both  phases  as  a  function  of  outer  iterations  (Fig.  22),  and  (ii) 
normalized  CPU  time  (Table  2)  needed  for  the  maximum  normalized  residual  of  all  variables  and  for 
all  phases  to  drop  below  £s“10'^.  In  order  not  to  overload  plots,  only  residuals  over  the  densest  grid 
using  the  SG,  PG,  and  FMG  methodologies  are  presented.  As  shown  (Fig.  22),  the  PG  and  FMG 
methods  greatly  reduce  the  number  of  iterations  with  the  number  required  by  the  PG  method  being 
around  four  times  the  number  required  by  he  FMG  method.  The  slopes  of  the  various  curves 
displayed  in  Fig.  22,  which  are  clearly  much  higher  for  the  FMG  method,  further  reveal  this. 

The  normalized  CPU-times,  required  by  the  four  algorithms  to  solve  the  problem  on  the  various  grids 
and  with  the  various  methodologies,  presented  in  Table  2  indicate  an  increase  with  increasing  the  grid 
size.  Moreover,  higher  savings  are  realized  with  both  the  PG  and  FMG  methods  over  the  ones 
obtained  in  the  previous  problem  over  all  grids.  Again  the  CPU-time  required  by  the  four  algorithms 
is  very  close  with  no  appreciable  superiority  of  any  one.  Using  the  PG  method  on  the  finest  mesh 
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accelerates  the  convergence  by  a  factor  of  3.6  whereas  the  FMG  does  a  much  better  job  of 
accelerating  the  solution  by  a  factor  of  around  9. 

Compressible  dilute  air-particle  subsonic  flow  over  a  flat  plate 

This  problem  has  been  used  by  several  researchers  [101-107]  to  check  their  computational 
methodologies  and  is  used  here  for  the  same  purpose.  It  is  known  that  two-phase  flow  greatly  changes 
the  main  features  of  the  boundary  layer  over  a  flat  plate.  Typically,  three  different  regions  are  defined 
in  the  two-phase  boundary  layer  (Fig.  23),  which  can  be  distinguished  by  the  relative  velocity 
between  the  two  phases:  a  large-slip  region  close  to  the  leading  edge,  a  moderate-slip  region  further 
down,  and  a  small-slip  one  far  downstream.  It  has  also  been  established  that  [104],  the  characteristic 
scale  in  this  two-phase  problem  is  the  relaxation  length  Xc,  defined  as: 


Eventhough  variations  in  gas  density  are  small  under  the  conditions  considered,  these  variations  are 
not  neglected  and  the  flow  is  treated  as  compressible  for  the  continuous  phase  and  as  incompressible 
for  the  dispersed  phase.  Moreover,  drag  is  the  only  interfacial  force  retained  due  to  its  dominance  over 
other  interfacial  forces.  Denoting  the  continuous  and  dispersed  phases  by  superscripts  (c)  and  (d). 


respectively,  this  force  is  computed  as  [104]: 


^  ’"p 


(133) 


(134) 
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where  the  drag  coefficient  is  given  hy: 

C„=— Re  +-Re"'’  (135) 

In  the  energy  equation,  heat  transfer  due  to  radiation  is  neglected  and  only  convective  heat  transfer 
around  an  isolated  particle  is  considered.  Moreover,  the  particles  have  no  individual  random  motion, 
mutual  collisions,  and  other  interactions  among  them.  Therefore,  only  the  process  of  drag  and  heat 
transfer  couple  the  particles  with  the  gas.  Under  such  conditions,  the  interfacial  terms  in  the  gas 
(continuous  phase)  and  particles  (dispersed  phase)  energy  equations  reduce  to  [104]: 


t(c)  _  1  p 

-^g-p+Vp-" 

(136) 

=  -0 

Vg-p 

(137) 

where: 

(138) 

1  1 

Nu=2.0+0.6Re2(Pr^"^)3 

(139) 

(140) 

In  the  above  equations,  Nu  is  the  Nusselt  number,  the  gas  Prandtl  number ,  the  gas  thermal 
conductivity,  T  the  temperature,  and  other  parameters  are  as  defined  earlier. 


In  the  simulation,  the  particle  diameter  is  chosen  to  be  10  pm,  the  particle  Reynolds  number  is 
assumed  to  be  equal  to  10,  and  the  material  density  is  1766  kg/m^.  The  Prandtl  number  is  equal  to 
0.75.  The  south  boundary  (wall)  is  treated  as  a  no-slip  wall  boundary  for  the  gas  phase:  both 
components  of  the  gas  velocity  are  set  to  zero,  while  the  particle  phase  encounters  slip  wall 
conditions.  The  normal  fluxes  are  set  to  zero.  The  gas  and  the  particles  enter  the  computational 
domain  under  thermal  and  dynamical  equilibrium  conditions.  A  mass  load  ratio  of  1  between  the 
particles  phase  and  the  gas  phase  is  used.  Results  are  displayed  using  the  following  dimensionless 
variables  in  order  to  bring  all  quantities  to  the  same  order  of  magnitudes: 

X*  =  y*  =  u*  =— ,  v*  =  — Re=^ — ^  (141) 

Figure  24  shows  the  results  for  the  steady  flow  obtained  on  a  rectangular  domain  with  a  mesh  of 
density  104x48  C.V.  stretched  in  the  y-direction.  The  figure  provides  the  development  of  gas  and 
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particles  velocity  profiles  within  the  three  regions  mentioned  earlier.  In  the  near  leading  edge  area 
(x*=0.1),  the  gas  velocity  is  adjusted  at  the  wall  to  obtain  the  no-slip  condition  as  for  the  case  of  a 
pure  gas  boundary  layer.  The  particles  have  no  time  to  adjust  to  the  local  gas  motion  and  there  is  a 
large  velocity  slip  between  the  phases.  In  the  transition  region  (x*=l),  sigmficant  changes  in  the  flow 
properties  take  place.  The  interaction  between  the  phases  cause  the  particles  to  slow  down  while  gas 
accelerate  as  apparent  in  the  plots.  In  the  far  downstream  region  (x*=5),  the  particles  have  enough 
time  to  adjust  to  the  state  of  the  gas  motion.  The  slip  is  very  small  and  the  solution  tends  to 
equilibrium.  These  results  are  in  excellent  agreement  with  numerical  solutions  reported  by  Thevand  et 
al.  [107]  plotted  in  Fig.  (24),  which  gives  credibility  to  the  proposed  methodology. 

Credibility  being  established,  the  problem  is  solved  over  four  different  grids  of  sizes  13x6,  26,12, 
52x24,  and  104x48  control  volumes  using  all  multi-fluid  algorithms  and  the  SG,  PG,  and  FMG 
solution  methods.  As  in  the  previous  two  problems,  results  are  displayed  in  the  form  of  total  mass 
residuals  (Fig,  25)  and  normalized  CPU  time  (Table  3)  with  8s=10’^  Only  residuals  over  the  densest 
grid  using  the  SG,  PG,  and  FMG  methodologies  are  presented. 

As  depicted  in  Fig.  25(a),  PISO  requires  the  least  number  of  iterations  for  the  reasons  stated  earlier 
(the  same  apply  for  PRIME  Fig.25  (f))  and  this  is  not  associated  with  the  lowest  computational  effort 
(Table  3).  Moreover,  the  convergence  behavior  of  SIMPLE  (Fig.  25(b)),  SIMPLEC  (Fig.  25(c)), 
SIMPLEM  (Fig.  25(d)),  and  SIMPLEX  (Fig.  25(g))  is  very  similar  requiring  nearly  the  same  number 
of  iterations  with  all  solution  methods.  In  comparison,  SIMPLEST  (Fig.  25(e))  consumes  a  larger 
number  of  iterations.  The  worst  performance  is  for  PRIME  (Fig.  25(1))  especially  with  the  SG  method 
where  it  is  seen  to  require  over  8000  iterations  to  reduce  the  total  mass  residuals  to  the  desired  level. 
The  above  can  also  be  inferred  from  the  plots  displayed  in  Fig.  25(h).  Again,  the  PG  and  FMG 
methods  greatly  reduce  the  number  of  iterations  with  the  number  of  iterations  required  by  the  PG 
method,  for  most  algorithms,  being  around  four  times  the  number  required  by  the  FMG  method. 

The  normalized  CPU-times  for  all  cases  considered  are  presented  in  Table  3.  The  general  trend  is 
similar  to  the  previous  cases  with  the  SG  requiring  the  highest  CPU  effort  and  the  FMG  being  the 
most  efficient.  However,  the  PG  method  is  not  as  efficient  as  in  the  previous  problem  whereas  the 
FMG  remains  reliably  competent.  In  fact,  the  use  of  the  PG  method  on  the  dense  grid  decreases,  on 


n  j  /  oiSJW/ c-A»ujcw  1  tiuic,  r  utwnK-  ivic,tnuu  vvun  r  utt.  ivjimn-\jt  lu  yroL-ctc;/ j%Jt  iviuin~i  lutu  i  iuvv  ut  Jtit  u//ccui> 


-r  / 


average  for  all  algorithms,  the  computational  time  by  about  37%  whereas  about  570%  reduction  in 
computational  time  is  realized  with  the  FMG  method.  Consistently,  the  CPU-time  required  by  the 
SIMPLE,  SIMPLEC,  and  SIMPLEX  algorithms  is  very  close  with  no  appreciable  superiority  of  any 
one  of  these  algorithms  over  the  others.  Moreover,  the  virtues  of  the  FMG  method  increase  as  the  grid 
size  increases  (e.g.  for  SIMPLEC  the  acceleration  rate  increases  from  1.2  with  two-levels,  to  2.26 
with  three-levels,  up  to  8.06  with  four  levels).  On  the  dense  grids,  which  is  of  interest,  the  most  time- 
consuming  algorithm  is  PRIME  followed  by  SIMPLEST  and  PISO.  The  best  performance  however,  is 
for  SIMPLE,  SIMPLEC,  and  SIMPLEX  with  that  of  SIMPLEM  being  close. 

Inviscid  transonic  dusty  flow  in  a  converging-diverging  nozzle 

As  a  final  test  for  the  newly  suggested  numerical  procedures,  dilute  two-phase  supersonic  flow  in  an 
axi-symmetric  converging-diverging  rocket  nozzle  is  considered.  Several  researchers  have  analyzed 
the  problem  and  data  is  available  for  comparison  [108-1 17].  In  most  of  the  reported  studies,  a  shorter 
diverging  section,  in  comparison  with  the  one  considered  here,  has  been  used  when  predicting  the 
two-phase  flow.  Two-phase  results  for  the  long  configuration  have  only  been  reported  by  Chang  et.al. 
[1 12].  The  flow  is  assumed  to  be  inviscid  and  the  single-phase  results  are  used  as  an  initial  guess  for 
solving  the  two-phase  problem.  The  physical  configuration  (Fig.  26)  is  the  one  described  in  [112]. 
The  viscosity  of  the  fluid  varies  with  the  temperature  according  to  Sutherland’s  law  for  air: 


=  1.458x10’ 


T^'^ -I- 110.4 


(142) 


The  coupling  between  gas  and  particle  phases  is  through  the  interfacial  momentum  and  energy  terms. 
The  force  exerted  on  a  single  particle  moving  through  a  gas  is  given  as  [1 13] 


f^-6jrr;,foli'''(v'‘''-v'‘=0  (144) 

so  that  for  N  particles  in  a  unit  volume  the  effective  drag  force  is 


where  fo  is  the  ratio  of  the  drag  coefficient  Cd  to 


(145) 

(146) 

the  stokes  drag  CDo=24/Rep  and  is  given  by  [1 12] 
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n,„7  0.0175Re„ 

fo  =  1  +  0. 15  Re“""+ - 4  ^  1,6 
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RCp  <  3x10 


The  heat  transferred  from  gas  to  particle  phase  per  unit  volume  is  given  as  [1 13] 

'X 

Qg_p  =  (148) 

2  tp 

Where  is  the  thermal  conductivity  of  the  gas  and  Nu,  the  Nusselt  number,  is  written  as  [1 13] 

Nu  =  2+0.459  RCp^  Prf  ^  (149) 

The  gas-particle  inter-phase  energy  term  is  given  by 

If  =  +-^ -  v^'>d  +|— X<'^Nu(t'''> (150) 

2  fp  2  Tp  2  Tp 

'X  r(d) 

if  =  (151) 

2  rp 

where  the  first  two  terms  on  the  right-hand  side  of  equation  (150)  represent  the  energy  exchange  due 
to  momentum  transfer. 

The  physical  quantities  employed  are  similar  to  those  used  in  [112].  The  gas  stagnation  temperature 
and  pressure  at  inlet  to  the  nozzle  are  555  °K  and  10.34x10^  N/m^,  respectively.  The  specific  heat  for 
the  gas  and  particles  are  1.07x10^  J/Kg®K  and  1.38x10^  J/Kg^IC,  respectively,  and  the  particle  density 
is  4004.62  kg/m^  With  a  zero  inflow  velocity  angle,  the  fluid  is  accelerated  from  subsonic  to 
supersonic  speed  in  the  nozzle.  The  inlet  velocity  and  temperature  of  the  particles  are  presumed  to  be 
the  same  as  those  of  the  gas  phase.  Results  for  two  particle  sizes  of  radii  1  and  10  |Ltm  with  the  same 
mass  fraction  (j)=0.3  are  presented  using  a  grid  of  size  188x80  C.V.  Figures  27(a)  and  27(b)  show  the 
particle  volume  fraction  contours  while  Figures  27(c)  and  27(d)  display  the  velocity  distribution.  For 
the  flow  with  particles  of  radius  1  |xm,  a  sharp  change  in  particle  density  is  obtained  near  the  upper 
wall  downstream  of  the  throat,  and  the  particle  density  decreases  to  a  small  value.  With  the  large 


particle  flow  (10  [im),  however,  a  much  larger  particle-free  zone  appears  due  to  the  inability  of  the 
heavier  particles  to  turn  around  the  throat  comer.  These  findings  are  in  excellent  agreement  with 
published  results  reported  by  Chang  et  al.  [112]  and  others  using  different  methodologies.  In  addition 
contours  are  similar  to  those  reported  by  Chang  et  al.  [112].  A  quantitative  comparison  of  current 
predictions  with  published  experimental  and  numerical  data  is  presented  in  Fig.  28  through  gas  Mach 
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number  distributions  along  the  wall  (Fig.  28(a))  and  centerline  (Fig.  28(b))  of  the  nozzle  for  the  one- 
phase  and  two-phase  flow  situations  with  particles  of  radii  10  pm.  As  can  be  seen,  the  one-phase 
predictions  fall  on  top  of  experimental  data  reported  in  [115-117].  Along  the  centerline  of  the  nozzle, 
current  predictions  are  of  quality  better  than  those  obtained  by  Chang  et  al.  [112].  Since  the  nozzle 
contour  has  a  rapid  contraction  followed  by  a  throat  with  a  small  radius  of  curvature,  the  flow  near  the 
throat  wall  is  overturned  and  inclined  to  the  downstream  wall.  A  weak  shock  is  thus  formed  to  turn 
the  flow  parallel  to  the  wall.  This  results  in  a  sudden  drop  in  the  Mach  number  value  and  as  depicted 
in  Fig.  28(b),  this  sudden  drop  is  correctly  envisaged  by  the  solution  algorithm  with  the  value  after  the 
shock  being  slightly  over  predicted. 

Due  to  the  unavailability  of  two-phase  flow  data,  predictions  are  compared  against  the  numerical 
results  reported  in  [112].  As  displayed  in  Figs.  28(a)  and  28(b),  both  solutions  are  in  good  agreement 
with  each  other  indicating  once  more  the  correctness  of  the  calculation  procedures.  The  lower  gas 
Mach  number  in  the  two-phase  flow  is  caused  by  the  heavier  particles  which  reduce  the 

gas  velocity.  Moreover,  owing  to  the  particle-tree  zone,  the  Mach  number  difference  between  the  one- 
and  two-phase  flows  along  the  wall  is  smaller  than  that  at  the  centerline. 

To  compare  the  relative  performance  of  the  multi-fluid  algorithms,  the  problem  is  solved  via  the  SG 
method  using  the  SIMPLE,  SIMPLEC,  SIMPLEX,  and  SIMPLEST  algorithms  over  three  different 
grids  of  sizes  47x20,  94x40,  and  188x80  C.V.  As  before,  results  are  displayed  in  the  form  of  total 
mass  residuals  (Fig.  29)  and  normalized  CPU  times  (Table  4)  with  Es=10‘^. 

As  shown  in  Fig.  29,  with  the  exception  of  SIMPLEST  on  the  94x40  grid,  all  algorithms  require 
almost  the  same  number  of  iterations  with  SIMPLE  and  SIMPLEC  being  slightly  better.  As  expected, 
the  number  of  iterations  increases  with  increasing  grid  density.  Moreover,  the  convergence  histories 
of  all  algorithms  are  nearly  identical.  The  normalized  CPU-times  presented  in  Table  4  confirm  these 
conclusions  and  reveal  the  close  performance  of  SIMPLE  and  SIMPLEC.  The  normalized  CPU-times 
needed  by  SIMPLEX  are  the  highest  due  to  the  additional  equations  solved.  Its  performance  could 
improve  on  denser  grids. 
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Closing  Remarks 

The  implementation  of  seven  MCBA  multi-fluid  algorithms  for  the  simulation  of  multi-fluid  flow  at 
all  speeds  was  accomplished.  The  algorithms  were  embedded  within  a  non-linear  full  multi-grid 
strategy.  A  two-fluid  k-e  model  and  several  inter-phase  models  were  also  employed.  Solving  a  variety 
of  one-  and  two-dimensional  two-phase  flow  problems  assessed  the  performance  and  accuracy  of 
these  algorithms.  For  each  test  problem,  solutions  were  generated  on  a  number  of  grid  systems  using 
the  single  grid  method  (SG),  the  prolongation  grid  method  (PG),  and  the  full  non-linear  multi-grid 
method  (FMG).  Results  obtained  demonstrated  the  capability  of  all  algorithms  to  deal  with  multi-fluid 
flow  situations  and  to  predict  multi-fluid  flow  at  all  speeds,  and  the  ability  of  the  FMG  method  to 
tackle  the  added  non-linearity  of  laminar  and  turbulent  multi-fluid  flows.  The  convergence  history 
plots  and  CPU-times  presented,  indicated  similar  performances  for  SIMPLE,  SIMPLEC,  and 
SIMPLEX.  The  PISO,  SIMPLEM,  and  SIMPLEST  algorithms  were  in  general  more  expensive  than 
SIMPLE.  In  general,  the  PRIME  algorithm  was  the  most  expensive  to  use.  The  PG  and  FMG  methods 
accelerated  the  convergence  rate  for  all  algorithms.  The  FMG  method  was  found  to  be  by  far  more 
efficient. 

It  is  the  author’s  view  that  the  accomplishments  are  immensely  greater  than  initially  anticipated 
however  there  remains  a  great  deal  of  work  yet  to  be  done.  The  items  demanding  additional  efforts  in 
the  immediate  future  are: 

1.  Implementation  of  the  GCBA  group  of  algorithms 

2.  Implementation  of  more  multi-fluid  models  (boiling,  condensation,  solidification,  . . .) 

3.  Implementation  of  additional  turbulence  models  (e.g.  k-co,  SST,. . .) 

4.  Implementation  of  multi-component  physics  for  the  individual  fluids 

5.  The  use  of  unstructured  grids 
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Table  1  Normalized  CPU-times  for  turbulent  bubbly  flow  in  a  pipe. 

Table  2  Normalized  CPU-times  for  turbulent  air-particle  flow  in  a  pipe. 

Table  3  Normalized  CPU-times  for  Dusty  flow  over  a  flat  plate. 

Table  3  Normalized  CPU-times  for  Dusty  flow  in  a  converging-diverging  nozzle. 
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Figure  Captions 

Fig.  1  (a)  Control  volume,  (b)  Typical  control  volume  faces  and  geometric  nomenclature. 

Fig.  2  Physical  domain  for  the  gas-particle  transport  problem. 

Fig.  3  (a)  Comparison  between  the  analytical  and  numerical  particle  velocity  distributions, 

(b)-(g)  convergence  histories  on  the  different  grid  systems,  (h)  and  convergence  histories  of 
the  various  algorithms  on  the  80  C.V.  grid  for  the  horizontal  dilute  gas-solid  flow  problem. 

Fig.  4  Normalized  CPU-times  for  the  horizontal  dilute  gas-solid  flow  problem. 

Fig.  5  (a)  gas  and  particle  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different  grid 
systems,  (h)  and  convergence  histories  of  the  various  algorithms  on  the  80  C.V.  grid  for  the 
horizontal  dense  gas-solid  flow  problem. 

Fig.  6  Normalized  CPU-times  for  the  horizontal  dense  gas-solid  flow  problem. 

Fig.  7  (a)  Liquid  and  gas  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different  grid 
systems,  (h)  and  convergence  histories  of  the  various  algorithms  on  the  80  C.V.  grid  for  the 
horizontal  dilute  bubbly  flow  problem. 

Fig.  8  Normalized  CPU-times  for  the  horizontal  dilute  bubbly  flow  problem. 

Fig.  9  (a)  Liquid  and  gas  velocity  distributions,  (b)-(g)  convergence  histories  on  the  different  grid 
systems,  (h)  and  convergence  histories  of  the  various  algorithms  on  the  80  C.V.  grid  for  the 
horizontal  dense  bubbly  flow  problem. 

Fig.  10  Normalized  CPU-times  for  the  horizontal  dense  bubbly  flow  problem. 

Fig.  1 1  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the  different 
grid  systems  for  the  vertical  dilute  gas-solid  flow  problem. 

Fig.  12  Normalized  CPU-times  for  the  vertical  dilute  gas-solid  flow  problem. 

Fig.  13  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the  different 
grid  systems  for  the  vertical  dense  gas-solid  flow  problem. 

Fig.  14  Normalized  CPU-times  for  the  vertical  dense  gas-solid  flow  problem. 

Fig.  15  (a)  gas  and  particle  velocity  distributions;  and  (b)-(h)  convergence  histories  on  the  different 
grid  systems  for  the  vertical  dilute  bubbly  flow  problem. 
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Fig.  16  Normalized  CPU-times  for  the  vertical  dilute  bubbly  flow  problem. 

Fig.  17  (a)  gas  and  particle  velocity  distributions,  and  (b)-(h)  convergence  histories  on  the  different 
grid  systems  for  the  vertical  dense  bubbly  flow  problem. 

Fig.  18  Normalized  CPU-times  for  the  vertical  dense  bubbly  flow  problem. 

Fig.  19  Comparison  of  fully  developed  liquid  velocity  and  void  fraction  profiles  for  turbulent  bubbly 
upward  bubbly  flow  in  a  pipe;  (a)  Seriwaza  et  al.  data,  (b)  Lahey  et  al.  data. 

Fig.  20  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid,  and 

(h)  convergence  histories  of  the  various  algorithms  on  the  finest  mesh  using  the  FMG  method 
for  turbulent  upward  bubbly  flow  in  a  pipe. 

Fig.  2 1  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  for  turbulent  air-particle 
flow  in  a  pipe. 

Fig.  22  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEST,  and  (d)  SIMPLEX 
algorithms  using  the  SG,  PG,  and  FMG  methods  on  the  finest  mesh  for  turbulent  air-particle 
flow  in  a  pipe. 

Fig.  23  The  three  different  regions  within  the  boundary  layer  of  dusty  flow  over  a  flat  plate. 

Fig.  24  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  inside  the  boundary  layer  at 
different  axial  locations  for  dilute  two-phase  flow  over  a  flat  plate 

Fig.  25  (a)-(g)  Convergence  histories  of  the  SG,  PG,  and  FMG  methods  on  the  finest  grid,  and 

(h)  convergence  histories  of  the  various  algorithms  on  the  finest  mesh  using  the  FMG  method 
for  dusty  gas  flow  over  a  flat  plate. 

Fig.  26  Physical  domain  for  the  dusty  gas  flow  in  a  converging-diverging  nozzle. 

Fig.  27  (a,b)  Volume  Fraction  contours  and  (c,d)  particle  velocity  vectors  for  dusty  gas  flow  in  a 
converging-diverging  nozzle. 

Fig.  28  Comparison  of  one-phase  and  two-phase  gas  Mach  number  distributions  along  the  (a)  wall 
and  (b)  centerline  of  the  dusty  flow  in  a  converging-diverging  nozzle  problem. 

Fig.  29  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEST,  and  (d)  SIMPLEX 
algorithms  using  the  SG  method  for  dusty  gas  flow  in  a  converging-diverging  nozzle. 


Table  1  Normalized  CPU-times  for  turbulent  bubbly  flow  in  a  pipe. 


Table  2  Normalized  CPU-times  for  turbulent  air-particle  flow  in  a  pipe. 


Table  4  Normalized  CPU-times  for  Dusty  flow  in  a  converging-diverging  nozzle. 


Fig.  2  Physical  domain  for  the  gas-particle  transport  problem. 


and  convergence  histories  of  the  various  algorithms  on  the  80  C.V.  grid  for  the  horizontal  dilute  gas-solid  flow  problem. 


Dilute  gas-solid  flow 


algorithms  on  the  80  C.V.  grid  for  the  horizontal  dense  gas-solid  flow  problem. 
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algorithms  on  the  80  C.V.  grid  for  the  horizontal  dilute  bubbly  flow  problem. 


Dilute  bubbly  flow 


Fig.  8  Normalized  CPU-times  for  the  horizontal  dilute  bubbly  flow  problem. 


Vertical  dilute  gas-solid  flow 


Fig.  12  Normalized  CPU-times  for  the  vertical  dilute  gas-solid  flow  problem. 
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Fig.  16  Normalized  CPU-times  for  the  vertical  dilute  bubbly  flow  problem. 
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Fig.  19  Comparison  of  fiilly  developed  liquid  velocity  and  void  fraction  profiles  for  turbulent 
upward  bubbly  flow  in  a  pipe;  (a)  Seriwaza  et  al.  data,  (b)  Lahey  et  al.  data. 
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finest  mesh  using  the  FMG  method  for  turbulent  upward  bubbly  flow  in  a  pipe. 


Fig.  21  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  for  turbulent 
air-particle  flow  in  a  pipe. 
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Fig.  22  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEST,  and 


(d)  SIMPLEX  algorithms  using  the  SG,  PG,  and  FMG  methods  on  the  finest  mesh  for 
turbulent  air-particle  flow  in  a  pipe. 


Region  III 


Fig.  23  The  three  different  regions  within  the  boundary  layer  of  dusty  flow  over  a  flat 
plate. 


Fig.  24  Comparison  of  fully  developed  gas  and  particle  velocity  profiles  inside  the  boundary  layer  at  different  axial  locations  for  dilute  two-phase  flow  over 


Fig.  26  Physical  domain  for  the  dusty  gas  flow  in  a  converging-diverging  nozzle. 
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Fig.  27  (a,b)  Volume  Fraction  contours  and  (c,d)  particle  velocity  vectors  for  dusty  gas 
flow  in  a  converging-diverging  nozzle. 
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Fig.  28  Comparison  of  one-phase  and  two-phase  gas  Mach  number  distributions  along  the  (a) 
wall  and  (b)  centerline  of  the  dusty  flow  in  a  converging-diverging  nozzle  problem. 
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Fig.  29  Convergence  histories  of  the  (a)  SIMPLE,  (b)  SIMPLEC,  (c)  SIMPLEST,  and 

(d)  SIMPLEX  algorithms  using  the  SG  method  for  dusty  gas  flow  in  a  converging- 
diverging  nozzle. 


